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ABSTRACT 
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Two theoretical methods and the development of a guidance, navigation, and 
control rapid protoyping system address the issue of considering the integral partici- 
pation of feedback early in the design process. The first method addresses the problem 
of sizing the horizontal tail on a statically unstable transport aircraft. Dynamic con- 
straints, which implicitly involve the flight control system, may prove more restrictive 
than traditional static constraints. Recovery from a severe angle of attack excursion, 
or penetration of a vertical wind-shear, are formulated in terms of the solution to 
a convex minimization problem utilizing LMIs, and numerically solved. The second 
method addresses the problem of tracking inertial trajectories, with applications for 
unmanned air vehicles. This problem is posed and solved within the framework of 
gain-scheduled control theory. This leads to a new technique for integrated guidance 
and control systems, with guaranteed performance and robustness properties. Finally, 
a rapid prototyping system for the flight testing of guidance, navigation, and control 
algorithms for unmanned air vehicles is designed. The system affords a small team the 
ability to take a new concept in guidance, navigation, and control from initial concep- 
tion to flight test. A proof-of-concept demonstration of the system is detailed when 
the new integrated guidance and control algorithm previously described is tested in 
flight on an unmanned air vehicle. 
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I. 



OVERVIEW 



A. LAYOUT OF THE CHAPTERS 

Methods in integrated plant-controller design and integrated guidance-controller 
design are contained in three chapters and two appendices. Applicable supporting 
computer code is explained in the appendices. The organization and content of the 
body of the work is as follows. 

Chapter II addresses the problem of sizing the horizontal tail on a statically 
unstable transport aircraft. A brief background in linear matrix inequalities prepares 
the reader for the problem formulation. Recovery from a severe angle-of-attack ex- 
cursion, or penetration of a vertical wind-shear, is formulated in terms of the solution 
to a convex minimization problem, and numerically solved. The effects of a num- 
ber of parameters in the aircraft definition, such as canards and flexible motion, are 
addressed. 

In Chapter III, the problem of unmanned air vehicles tracking inertial trajec- 
tories is addressed. The concept of trimming trajectories is defined, and the tracking 
of a trimming trajectory is converted into the problem of driving a generalized er- 
ror vector - which implicitly includes the distance to the trajectory - to zero. The 
linearization of the generalized error dynamics along the trajectory is shown to be 
time-invariant. Using these results, the problem of trajectory tracking is posed and 
solved within the framework of gain-scheduled control theory. This leads to a new 
technique for integrated guidance and control system design for unmanned air vehi- 
cles. 

The development of a rapid prototyping system for flight testing of guidance, 
navigation, and control algorithms for unmanned air vehicles is presented in Chapter 
IV. The system affords a small team the ability to take a new concept in guidance, 
navigation, and control from simulation to flight test. A proof-of-concept flight test 
demonstration of a new integrated guidance and control algorithm is detailed. The 
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project includes a parameter identification problem for an unmanned air vehicle at the 
Naval Postgraduate School. Controller synthesis and implementation follow, based 
on the identified model. Finally, the system is used to test the theory presented in 
Chapter III in flight. 



B. CONTRIBUTION 

This section states the original contribution of this work. 

• The problem of sizing a horizontal control surface for a statically unstable 
transport aircraft is addressed. Closed-loop performance of the plant and 
controller in response to certain dynamic constraints is used to define a design 
region that captures the integrated nature of the plant definition and controller 
synthesis process. 

• A numerical tool is developed that determines acceptable combinations of 
tail volume, center-of-gravity location, and peak actuator rate for a statically 
unstable transport aircraft using state-feedback, or output-feedback, linear 
control to recover from a severe angle-of-attack excursion, or vertical wind- 
shear penetration. The method formulates the dynamic constraint as a convex 
optimization problem, and uses recently developed LMI solvers to solve the 
problem numerically. 

• Numerical results demonstrate how to quantify the effect of the following in 
terms of their influence on the sizing of the horizontal control surfaces. 

- The addition of a canard to an aircraft definition. 

— The use of dynamic output versus static, state-feedback linear control. 

— The effect of simple, symmetric aeroelastic motion. 

• Work on implementing nonlinear gain-scheduled controllers has been extended 
to trajectory tracking control problems for unmanned air vehicles. A broad 
class of trimming trajectories is defined, and a technique is developed that 
provides independent control of the space and time coordinates along trajec- 
tories in the class. Using this technique, it is shown that the performance and 
robustness of the linear design is recovered locally by the nonlinear plant and 
gain-scheduled controller. 

• A unified system for the synthesis, simulation, and flight testing of avionics 
algorithms on unmanned air vehicles is constructed. The system is used to 
take new theory in integrated guidance and control from conception to flight 
test. 
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II. INTEGRATED PLANT CONTROLLER 

OPTIMIZATION 

A. INTRODUCTION 

The use of an automatic flight augmentation system is commonplace on a 
modern aircraft. The benefits of its use may include such things as the remedy of 
undesirable flight characteristics, the reduction of pilot workload, and an increase in 
performance and fuel efficiency. Whether required for safe flight or not, it is hard 
to imagine a new transport aircraft being built without a sophisticated flight aug- 
mentation system. In fact, the next generation of high-speed civil transport (HSCT) 
aircraft will require some form of feedback control for safety of flight considerations, 
since current designs have an unstable short period. 

Since HSCT will require some form of automatic flight control always being 
active, the sizing of its control surfaces is no simple matter. Traditionally, static con- 
straints have been used to size the horizontal tail. For a given tail volume, constraints 
are calculated that limit the fore and aft travel of the center of gravity. Constraints 
that limit the forward center-of-gravity position include (1) sufficient nose-up pitch 
acceleration at the rotation speed (nose- wheel lift off), and (2) sufficient nose-up 
pitch acceleration at the approach speed in the landing configuration (go- around). 
Constraints that limit aft center-of-gravity position include (1) at brake release with 
maximum thrust, sufficient weight on the nose gear (tip back), (2) pitch-up acceler- 
ation at the rotation speed (nose-wheel lift off), and (3) sufficient nose-down pitch 
acceleration at minimum flying speeds [Ref. 26]. 

At the aft center-of-gravity locations projected for the approach flight condi- 
tion of the HSCT, dynamic constraints may be more restrictive than static constraints. 
The need to include dynamic considerations in the configuration definition process 
has been addressed before. References [Ref. 5] and [Ref. 4] describe early published 
work by Beaufrere in this area, largely motivated by the X-29 research program. In 
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[Ref. 41], Schmidt uses the system sensitivity function to describe the fundamental 
trade-off that exists between the level of static instability that can be controlled and 
vehicle flexibility. Previous (unpublished) work in industry has used time domain 
analysis to determine how far aft (how unstable) the center of gravity could be before 
the airplane was unrecoverable. The analysis included a given design angle-of-attack 
disturbance, and given rate and position limits of the pitch control effector. 

This work extends the work of [Ref. 26], and uses a numerical technique similar 
to previous work in [Ref. 33]. There an integrated aircraft controller design method- 
ology using LMIs was applied to the control power sizing for an F-14 aircraft. The 
major contribution of this work is two fold. First, the tail sizing design problem is 
defined in terms that include the integral participation of a feedback control system. 
Since the degree of control of the longitudinal dynamics depends on how fast and far 
the longitudinal control surface(s) can be moved by the control actuator(s), a natural 
metric that captures the size of the automatic flight control system is the maximum 
actuator rate. The design trade-off naturally includes consideration of actuator per- 
formance. For instance, it may be more cost effective to incorporate faster, generally 
larger and more expensive, actuators rather than pay the drag penalty associated 
with a larger horizontal tail. In that light, we introduce the following definition. 



Tail Sizing Design Space: 

The Tail Sizing Design Space is the region of “acceptable” combinations of tail 
volume ( Vh ), center- of- gravity station ( x c , g ), and peak actuator rate ( u max ). 
The triplet { Vf/, x c . 5 ., u max }, defines an aircraft model and an automatic flight 
control system . The model is obtained through the linearization of the nonlin- 
ear dynamics of the HSCT at an equilibrium point . It is partially defined by 
the specified tail volume and center- of- gravity position. The automatic flight 
control system is characterized by the specified maximum actuator rate in the 
triplet . By “ acceptable ” it is meant that , for the model associated with the 
triplet {V//, x c#5> , u m ax}; a linear controller is known to exist that 1) stabilizes 
the plant, 2) meets prescribed dynamic performance constraints, a7id 3) does 
not exceed the maximum actuator rate in response to the dynamic constraint . 
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Second, a numerical tool is developed that determines the Tail Sizing Design 
Space for a given HSCT dynamic model, flight condition, and dynamic constraint. 
This tool is termed the Tail Sizing Design 7oo/, and provides the capability to measure 
the effect of adding a second horizontal control surface in the form of a canard. The 
Tail Sizing Design Tool also provides the capability to measure the effect of simple, 
symmetric, flexible motion of the vehicle. Based on numerical analysis of designs, 
conclusions are drawn as to the relative effective of the use of canards, the use of 
static versus dynamic controllers, and the inclusion of aeroelastic effects. 

1. Example of a Design Space 

As a motivating example, consider a simple case. Let a(s) and b(h) be two in- 
dependent variables that are smooth functions of s and /i, respectively. Furthermore, 
assume that they are independent of time, i.e. — 0 and = 0. Let 

dx 

— — a(s)x(t) + b(h)u(t), (II. 1) 

at 

where x, u 6 7Z 1 . Then, equation II. 1 describes a family of dynamic systems. 

Evaluating (II. 1) about arbitrary values of the parameters a(s) and b(h ) results in 
the LTI systems 

x(t) = a(s 0 )x(t) + b(h 0 )u(t), 

x(0) = x 0 . (II. 2) 

It is desired to use feedback of the form 

u(t ) = — k(so, ho)x(t) (II. 3) 

to stabilize (II. 2). Lyapunov stability theory states that this is equivalent to the 
existence of a p > 0 such that 

2p{a(s 0 ) — b(h 0 )k(s 0 ,h 0 )} <0. (II.4) 

This implies that 

> km <IL5) 
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guarantees local stability of (II. 1) at (s 0 , h 0 ). This places a lower bound on k(s 0 , h 0 ). 

Assume that the absolute value of u(t) is constrained to be less than some 
nominal value, u max . For this simple system, the maximum value of |u(t)| occurs at 
t = 0 and is equal to 

Umax — j /i. (.Sq , M*o|- (H.6j 

This places an upper bound on k(s 0 , /io), 

k(s 0 ,h 0 ) < — . (II. 7) 

Xq 

The two limits, (II. 5) and (II. 6), are shown graphically in Figure 1. 



Region of 
Acceptable Gains 






g (A 

m 



u(max) 

xO 



Figure 1 . Range of Feedback Gains 



k(s, h) 



If the boundaries, are allowed to approach each other, the region 

of acceptable feedback gains approaches a single point where 

a(s 0 ) 



Ur 



b(h 0 ) 



x 0 



(II-S) 



Suppose a(s) varies linearly from a minimal value of 0.1 to a maximum value of 1.0. 
Suppose b(h) varies linearly from a minimal value of 0.3 to a maximum value of 0.9. A 
value of 1 was used for x 0 ■ The family of LTI systems, which validate the relationship 
in expression II. 8, is represented by the curved surface shown in Figure 2. 

Below the surface, the LTI systems are not stabilizable subject to the stated 
constraint. Above the surface, multiple controllers exist that stabilize a given system 
subject to the stated constraint. If the maximum acceptable value of u(t) is included 
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3.5-1 




Figure 2. Family of LTI systems with a single acceptable gain. 



as a horizontal plane, the Design Space for this simple triplet, {a, ,6 ,u maj J, consists 
of the region above the surface in Figure 2 and below this horizontal plane. This is 
shown in Figure 3. 




Figure 3. Example of a Design Space 
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The intersection of the two surfaces is a line which represents limiting combi- 
nations of a(s) and b(s) for a fixed u max and x 0 . 

The chapter begins with a background in Linear Matrix Inequalities sufficient 
to formulate the general problem in a form suitable for solution numerically. Next, 
the Tail Sizing Design Tool is developed for a rigid-body HSCT model. The dynamic 
constraint associated with the problem is a severe angle-of-attack disturbance. Nu- 
merical results follow based on an application of the Tail Sizing Design Tool to an 
HSCT model representative of current designs. Then, a method is developed that 
models the effects of simple, symmetric, flexible motion of the HSCT. Numerical re- 
sults follow based on applying the Tail Sizing Design Tool to an aeroelastic HSCT 
model. Comparisons to the rigid-body model are made. Finally, the Tail Sizing De- 
sign Tool is developed for a rigid-body HSCT model where the dynamic constraint 
is the penetration of a vertical wind-shear. The chapter ends with conclusions and 
recommendations. 

B, BACKGROUND IN LINEAR MATRIX INEQUALI- 
TIES 

It is usually not possible to express the solution to a complex problem with 
multiple constraints analytically. Progress in computational power, however, make 
numerical solutions to these problems attractive. Efficient convex optimization algo- 
rithms, recently made available from The Math Works [Ref. 17], allow problems that 
can be formulated in terms of LMIs to be solved exactly. 

Two major results in LMI theory are drawn on to formulate the problem of 
sizing a longitudinal control surface. The first concept is that of an invariant ellipsoid. 

Consider an LTI system: 



x — Fx , 

* = Gx, (II. 9) 
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where x G 7Z n , z £ 7Z P . Let P > 0 and define 

£ = {c e n n : ( T P( < 1}. 

Then £ is called an invariant ellipsoid associated with (11.9) if, for every trajectory x 
satisfying (II. 9), a’(O) € £ implies x(t) is in £, for alH > 0 [Ref. 6]. 

Let £ be an invariant ellipsoid for (II. 9), and define R(x) = x T Px. From the 
definition of £, it follows that V(x) < R(;r(0)) < 1, and 

V(x) = x T (F T P + PF)x < 0, 

=► F T P + PF < 0. 

On the other hand, suppose 3P > 0 such that F T P + PF < 0 and x(0) T Px(0) < 1. 
Then, for x(t) satisfying (II. 9), we obtain 

x t (F t P + PF)x =: V{x) < 0. 

Integrating both sides from 0 to T > 0 we get 

V(x(T))-V(x(0)) < 0, 

^ V(x(T)) < K(ar(0)) < 1, 

for any T > 0. Therefore, we have shown that £ is an invariant ellipsoid associated 
with the linear system (II. 9) if and only if F T P + PF < 0. 

Next, we will show how the idea of invariant ellipsoids can be. used to obtain 
bounds on the peak of the output in response to a known initial condition. Let £ be 
an invariant ellipsoid associated with the linear system (II. 9). Then, 

||z(£)|| 2 = z(t) T z{t) < ma x( T G T G( Vt > 0. 

Now, 

ma x( t G t G( = ma x u T P~ 1/2 G T GP~ 1/2 u, 

Cee IMI<i 

= II GP-^IW 

= A max (P- l/2 G T GP- 1 ' 2 ), (11.10) 
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where A max (.) denotes the maximum eigenvalue of the argument. Therefore, an upper 
bound on the peak of ||z(/)|| i n response to x(0) can be found as a square root of the 
minimum of 8 subject to: 



8- \ max (P- l ' 2 G T GP- 1/2 ) > 0, 


(mi) 


x(0) T Px(0) < 1, 


(11.12) 


F t P + PF < 0. 


(11.13) 



Nonlinear inequalities, such as equation 11.27, can be converted to LMI form using 
Schur complements [Ref. 22]. In general, any nonlinear inequality of the form 

R> 0, Q-SR~ 1 S T > 0, (11.14) 



for symmetric Q and R, is equivalent to the LMI 



Q s 

S T R 



> 0 . 



(11.15) 



To see this, consider the congruence transformation, 



I -SR- 1 




’ Q 


S 




I 0 




Q - SR' 1 S T 


0 


0 / 




_S T 


R 




-R~ T S T I 




0 


R 



, (11.16) 



where the block diagonal structure makes equivalence to expression 11.14 obvious. 
Using Schur complements on expression 11.27 yields, 



8-\ max (P- 1 ' 2 G T GP- 1 ' 2 ) > 0 «=* 



81 - P~ 1/2 G T GP~ 1/ 2 > 0 4=» 



8 GP- 1 ' 2 




I 0 




8 G 




/ 0 


P- 1 ' 2 G T / 




a, 

o 

1 




i 

s 

i 




0 P - 1 ' 2 



8 


G " 


>0 4=^ 


P 


G T ' 


G T 


P 




G 


8 



In summary, the square root of a solution to the following optimization problem: 
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min 8 



subject to 

P G T 

> 0, (11.17) 

G 8 

x(0) T /Lr(0) < 1, (11.18) 

F t P+PF < 0, (11.19) 

P > 0, (11.20) 



provides an upper bound on the peak of ||z(£)|| f° r the system (II. 9) in response to 
the initial condition x(0). 

The second concept used is that of an LMI region. In their recent work, 
Gahinet and Chilali [Ref. 7] introduce the concept of an LMI region, and have shown 
that such regions can be characterized by LMIs. LMI regions include sectors, disks, 
conic sectors, strips, as well as the intersection of any of these. For our purposes, the 
intersection of just two regions will suffice. The first region is a conic sector centered 
at the origin with half inner angle <f>. The second region is the left half-plane with its 
right boundary at — /3. The intersection of these regions is shown as the region C in 
Figure 4. 




Re(s) 
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The eigenvalues of the system described by equation II. 9 lie in the region C , 
if and only if there exists a positive definite matrix P, such that 



sin <f>(F T P + PF) cos <f>(-PF+F T P) 0 

cos (f>(PF-F T P) sin (f>{F T P+PF) 0 

0 0 F T P + PF + P2p 



< 0 . ( 11 . 21 ) 



1. State Feedback Synthesis LMIs 

A feedback interconnection of an LTI system and state-feedback controller is 



x = Ax + Bu 
* z = Cx + Du 
u = K x 



x — [A + BI\]x 
z = [C + DK]x, 



( 11 . 22 ) 



where x 6 lZ n , u G 1Z m , z € W . Substitution into the previous invariant ellipsoid 
and pole placement region matrix inequalities result in: 



111 1 0 Zrriax 

subject to 



P (C + DK) T 
(C + DK) 4„ 

x(0) T Px(0) 



(sin <t>)(A T P + K T B T P+ PA + PBK) (cos 4>){-PA - PBK + A T P + K T B T P) 

(cos <t>)(PA + PBK - A T P - K T B T P) (sin 4>){A T P + K T B T P + PA + PB K) ... 
0 0 

0 

0 

A T P + K T B T P A PA A PBK + P2(3 



> 0, 
< 0, 

< 0 , 



P > 0. (11.23) 



The above matrix inequalities are no longer affine in the unknown parameters, 
since the state-feedback controller, K, and positive definite matrix, P , appear as 
multiplicative terms. Since P is positive definite, let Y = P~ A . Perform the following 
congruence transformations, and perform the substitution of variables, IT = I\Y . 
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Y 0 




0 I 





P (C + DI<) T 
( C + DI<) 

P (C + DI\) T 



„2 

^ max 



Y 0 

0 I 

Y ((C + DK )Y) T 
(C + DK)Y zi ax 

Y ( CY + DW) T 

(CY+DW) z* max 



0 ^ 



>0 4* 



0 4* 



> 0. 



(11.24) 



Y 0 0 
0 Y 0 
0 0 Y 



(sin 4>)(A T P + K T B T P + PA + PBK) (cos <t>){-PA - PBK -f A T P + K T B T P) 

(cos <p)(PA -f PBK — A T P — K T B T P) (sin <j>){A T P + K T B T P + PA+ PBK) ... 
0 0 

0 
0 

A T P + 1< T B T P + PA + PBK + P2/3 

(sin 4>)(Y A T + W T B T + AY + BW) (cos <*>)( -AY - BW + YA T + W T B T ) 

(cos <fi)(AY + BW - YA T - W T B T ) (sin <j>)(Y A T + W T B T + + BW) ... 

0 0 

0 
0 

YA T + W T B T + AY+ 5^ + Y"2/3 



< 0. 



y o 
o r 

o o 



0 



(11.25) 



Using Schur complements: 



x(0) T Px(0) < 



a:(0) T y- 1 a:(0) 

1 x(0) T 

a:(0) Y 



< 1 
> 0. 



(11.26) 



Once again, the expressions 11.24, 11.25, 11.26 constitute an LMI in the un- 
known matrix variables, W and Y . Suppose 3 Y > 0 and W such that they validate 
expressions 11.24, 11.25, 11.26. Then, K = WY~ l is the state-feedback controller that 
maintains ||z|| < z max in response to the initial condition x(0),Vt > 0. 
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2. Output Feedback Synthesis LMIs 

A feedback interconnection of an LTI system and strictly proper, output- 
feedback controller is 



x = Ax + Bu 
y = Cx 

< z = C z x + D z u 




Xfc A^xfo T Bky 



A BC k 
B k C A k 



V 



C z D z Dk 



(11.27) 



U — CfcXfc 

where x € 7 Z n ,x k 6 u G 1l m , z G 7 l p , y € 1l j , tj G 7 l 2n , and 
where y 1 = [x T xj]. The initial condition is r/(0) T = rfi = [x^ x]. o ]. 

The four analysis matrix inequalities based on using the concept of an invariant 
ellipsoid and pole placement region are gathered for convenience. 



P G T 
G 81 



> 0, 



P > 0, 

XqPx o < 1, 



(11.28) 

(11.29) 

(11.30) 



(sin <f>)(F T P + PF) 
(cos <j>)(PF- F T P) 

0 



(cos <f>)(—PF + F T P) 
(sin 4>){F t P + PF) 

0 



0 

0 

F T P + PF + P2(3 



< 0. 
(11.31) 



If the closed-loop system matrices in expression 11.27 were substituted into 
expressions 11.28, 11.29, 11.30, and 11.31, the expressions would not be affine in the 
unknown matrix variables. Similar to the state-feedback case, the key is to find a 
congruence transformation with a “linearizing” effect on the unknown matrix vari- 
ables. In a recent work by Scherer and Gahinet [Ref. 8], the authors prove that such 
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a transformation exists for the output feedback case. Consider the following partition 
of P and P -1 , 



P 



p-i 



S N 
N T Q 

R M 
M T V 



( 11 . 32 ) 

( 11 . 33 ) 



Since P and P 1 are both symmertic, so are R and S. 



From the identity 



we can infer 



or 



Also, 



and 



P~ l P = 



I 0 
0 / 



R M 
M t F 



5 N 
N T Q 

RS + MN T RN + MQ 
M T S + VN T M T N -f VQ 

RS + MN t = /, 

NM t = I -SR. 



1 


s 




I 




n t 




0 



R 




I 


M T 




0 



( 11 . 34 ) 



( 11 . 35 ) 



( 11 . 36 ) 



( 11 . 37 ) 
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These identities are used to define the key matrices used in the congruence transfor- 
mation as follows: 



and 



Notice that 



/ 




/ 5 


1 

o 




O 

s 

1 



=: X 2y 



1 


I S 




R I 




i 

o 

s 




1 

O 

E-i 

i 



=: X x . 



AVV = 



CO 

i 




r 

co 


0 N T 




i 

o 

> 

s 

* 



-i -i 



P, 



= P. 



Using X\ and A^, consider the following congruence transformation: 



(11.38) 



(11.39) 



(11.40) 



i 

o 

X 

i 




P G T 




X 1 0 




Xj PX x X ?G T 


0 I 




G 81 




o 7 




i 

0 
* 

c* 

1 



Expanding the (1,1) block of expression 11.41, we obtain 



(11.41) 



A'.'P.Y, = XfX,X{'X u 

— Y 
— A 1 A 2 , 



R M 




1 S 


I 0 




O 

> 

H 

1 



/ 

/ 5 



(11.42) 
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Expanding the (2,1) block of expression 11.41 for the general case of output feedback, 
vve obtain 



GX x 



C z D z C k 



R 1 
M t 0 



C z R + D z C k M T C z • 



The matrix inequality in 11.28 becomes 



R I MCI D J + R T CJ 

I S CJ 

C z R + D z C k M T C z 4a* 



> 0 . 



After a congruence transformation with A"i , the LMI in (11.29) becomes 



R I 
1 S 



> o. 



This LMI is accounted for as a principle sub-matrix of the LMI in (11.44). 
complements are used to express the LMI in (11.30) as 



i 4 

VO P- X 



> 0 . 



Assume x ko = 0, and consider the following congruence transformation: 



1 0 




1 


So 




1 0 




1 

£ 

1 


1 

o 

* 
to ^ 




. vo 


p-i 




1 

* 

o 




Xfvo Xjp-'X 2 



1 Vo -^ 2 
Xj Vo XJX, 



(11.43) 



(11.44) 



(11.45) 

Schur 



(11.46) 



(II.4T) 
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Expanding the (1,2) block of equation 11.47 yields, 



T 

Vo 



I s 

0 N T 



x T o XqS 



(11.48) 



After the appropriate congruence transformation, and substituting (11.48) into (11.47), 
we have 

1 Xq Xq S 

x 0 R I >0. (11.49) 

S T x 0 I S 

To see how the pole placement region LMI is transformed, look at the trans- 
formation of the (1,1) block of equation 11.31. 



F t P + PF < 0 



P~ 1 F T + FP~ l < 0 • 4 = 4 - 

xTp~ 1 f t x 2 + xTfp~ 1 x 2 < 0. 



Since P 1 = X\ X 2 1 and P T = X 2 T Xl [ , equation 11.50 becomes 



(11.50) 



X?F T X 2 + XjFX, <0. 
For output feedback synthesis 



(11.51) 



F = 



A BCk 
B k C A k 

Substituting (11.52) into the first term in (11.51) results in 

A BC k 



Xl FAT = 



I 


0 




s 


N 





B k C Ak 





R I 




i 

O 

C 

i 



(11.52) 
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1 

Q 

?r 

I 




R / 


SA + NB k C SBC k + NA k 




M T 0 



AR+BC k M T A 

SAR + NB k CR + SBC k M T + NA k M T SA + NB k C 

(11.53) 

Define the change of variables, 



c k = 


C k M T , 


(11.54) 


B k = 


NB k C, 


(11.55) 


A k = 


SAR + NB k CR + SBC k M T + N A k M T , 


(11.56) 



and equation 11.53 is rewritten as 



AR+BCk A 
A k SA + B k 



(11.57) 



Each term in the pole placement LMI is transformed in a similar fashion. This 
results in 



sin <^(— T + E) cos <f>(— p + E T ) 
cos <f>(E — E T ) sin <j)(E T + E) 



0 

0 



0 



0 



l T + E + 



< 0 . 



R I 
I S 



(11.58) 



2/3 



Notice that expressions 11.49, 11.45, 11.44, and 11.58 are affine in the unknown 
matrix variables A k , B k , C k , R, and S. Thus, they constitute an LMI. Given 
A k , B k , C k , R, and S that validate expressions 11.49, 11.45, 11.44, and 11.58, the 
output-feedback controller is reconstructed as follows. First, reconstruct the matrices 
M and N via a singular value decomposition. 
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UI?V T 



I - RS, 

sv, 

i/s. 



(11.59) 

(11.60) 
(11.61) 



TV " 1 = 

M~ t = 

Form the controller by inverting the change of variables in expressions 11.55, 11.56, 
and 11.56 as 



A k 


= N~ 1 (A k - SA P R - B k C p R - SB p C k )M ~ T , 


(11.62) 


B k 


= N~'B k , 


(11.63) 


c k 


= C k M' T . 


(11.64) 



C. HIGH ANGLE OF ATTACK RECOVERY 

1. Problem Formulation for a Rigid Body HSCT 

The Plant Controller Optimization (PCO) problem to be solved in this section 
can be stated as follows: 

Let an HSCT flight condition be specified by the aircraft’s flight speed, altitude , 
and flight path angle . Furthermore , let a certain set of flying quality require- 
ments exist for the HSCT at that flight condition. Let the dynamic constraint 
be defined as the recovery of the aircraft from a high angle- of- attack excursion, 
while not exceeding certain peak actuator rate and actuator amplitude limits . 
Then, for a range of horizontal tail volumes, and for a range of peak actua- 
tor rates, determine the aft center- of- gravity limits for which there still exists 
a linear feedback controller that 1) stabilizes the plant, 2) satisfies the flying 
quality requirements, 3) satisfies the dynamic constraint. 

Let x cg denote the center-of-gravity location as a fraction of the reference 
chord, and let x denote the vector of HSCT longitudinal states. Let u denote the 
horizontal tail incidence angle, and let Vh denote the tail volume. Figure 5 provides 
a description of the tail volume, Vh- 

This representation is convenient since, for one flight condition, the wing-body 
mean aerodynamic center remains constant. Thus, the control surface volume and 
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center of gravity are decoupled in terms of their influence on the vehicle’s dynamics. 
The longitudinal dynamics of the aircraft can be expressed by the following system 
of differential equations: 



{ X = F(x,u,x cg ,V H ), 
z = H(x, u), 



(11.65) 



where JF(.) is a C 1 function of x,u,x cg , Vh, 7i is a C 1 function relating x and z, and 
z = [Vj 7 ] T defines the true airspeed ( V ? ) and flight path angle ( 7 ). Let x 0 ,u 0 denote 
the trim values of x and u for a given z 0 and center-of-gravity location, x cgo , i.e. 
T(xq ,u 0 ,x cgo , Vjf 0 ) = 0 and H(xo,u 0 ) = z 0 . Here the vector z is used to characterize 
the flight condition. Later, when x cg is allowed to change, the trim values, Xo and uo, 
will be recomputed for given values of x cgo , z 0 and V// 0 . Now, linearizing Q about 
Xq, uq yields a system of linear differential equations, 



Sx v4.(x.o, %cgo , f Hq )<^X “I - /?(Zq, X cgQ , V// Q )<5u, 



( 11 . 66 ) 
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where 8x and 8u denote small perturbations in x and u about xo and ti 0 , respectively. 
In the numerical analysis, the velocity components of 8x were normalized by the trim 
value of the true airspeed. 

Let the actuator dynamics be independent of (Zo, x cso , V// 0 ), and let them be 
given by the set of differential equations, 



8 x d d d 8 x d I ^3 1 d 8 w , 

8u = C a 8x a , (11.67) 

8u = C r 8x a + D r 8u, 

where 8u denotes the actuator amplitude and 8u denotes the actuator rate. Append- 
ing the actuator dynamics in series with the linearized longitudinal dynamics results 
in the system: 



Ql ("0? %cgo i f Hq ) ( 



8p = 

8u = 
8u = 



A BC a 
0 A a 
0 C a 
0 C T 



8v + 



0 

B a 



8u 



8u 

8v + D r 8u , ' 



( 11 . 68 ) 



where 



8u t = 



8x T 8xT 



(11.69) 



Flying quality requirements are typically characterized by the level of attention 
and skill required of the pilot to control the aircraft. They are grouped in three levels. 
A lower level corresponds to more benign flight characteristics. In order to achieve 
certain Level II flying qualities requirements, the eigenvalues of Q\ must be placed 
in a more restrictive region in the left half plane. Figure 6 shows suggested locations 
of the Category B, closed-loop pole locations of the HSCT short-perod mode. These 
locations meet Level II flying quality requirements for the flight condition used in this 
study. Notice, these locations can be characterized as having a minimum damping 
ratio of 0.2 and natural frequency of 0.2 radians per second. This suggests that, in 
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order to meet Level II requirements for the HSCT, the short-period eigenvalues of 
the closed-loop system Qi must be placed in the region in the complex plane given 
by Figure 7. 




Figure 6. Acceptable Short Period Pole Locations 
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2. Canard Configuration 

The horizontal tail is the most significant control effector for longitudinal con- 
trol of the aircraft. It is not uncommon, however, for long, slender aircraft designs 
to include canards for supplemental longitudinal control. The XB-70 and B-1B are 
two examples. Their presence is usually attributed to flying quality control issues 
involving the flexible nature of these aircraft [Ref. 1]. This point is addressed in a 
subsequent section concerning control of an aeroelastic aerodynamic model. To lay 
the ground work for that section, and as a baseline point of comparison for a rigid 
structure, this section addresses the addition of a longitudinal control effector for- 
ward of the wing. The effect on the Tail Sizing Design Space is explored when the 
feedback control system is free to utilize both control surfaces, in order to recover 
from angle-of-attack excursions. 

The revisions to the nonlinear equations of motion required to account for 
the addition of the canard parallel the development in [Ref. 14]. When these non- 
linear equations of motion are linearized about the equilibrium point determined by 
(zfhXcgo, V// 0 > V(7o)> the following LTI system results, 



8x — Xcgo ? ^ Ho ? Vco T D c *^cgo ? ^Co T (^ 0 ? *^cgo ? ^ Ho )^^/i ? (11.70) 

where u c and Uh represent movement of the canard and horizontal tail, respectively, 
and Vcq is the canard volume. Actuator dynamics for both the canard and horizontal 
tail are appended to the linearized longitudinal dynamics as before. The actuator 
dynamics are given by: 
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Aa c 


0 




Ba c 


0 




8u c 


8x a 


— 






bx a + 














0 


Aa h 




0 


B ah 
















T 






. _ 


8x a 


— 




5 s 











8u c = CaJXa, 

Suh C a ^8 X (j , 

8u c — C Tc 8x a -f- D Tc 8u c , 
8u h — C^ rh 8x a -I - D Th 8xih- 

This results in the following system, 



Ql (^0? %cgo ? ^ Ho ) ^Co ) ^ 



88 = 

8u c = 
8u c — 
8uh = 

8u h = 



where 



8v = 



A B c C ac BhC ah 
0 A ac 0 

0 0 A a 

0 C ac 0 



*<»/> 



8v + 



0 C Tc 0 



o o c, 



o-h 



o o a 



Th 



8v 

8v + 
8v 
8v -f 



D Tc 0 



l T 



8x T 8xJ 8x T" 

a c o>h 



B ae 0 




8u c 


0 B ah 




8u h 



8u c 

8uh 



r -i 




8u c 


0 D Th 




8u h 







(11.72) 



The initial condition is defined by the angle-of-attack perturbation. Let the 
angle-of-attack perturbation be given as ao; then the initial condition is 

■ \T 



8v o — 



0 V t cos 1 (ao) 0 0 0 0 



For each design point, Qi ( zo,x cgo , V// 0 , Vc 0 ), the question becomes, is the set of feed- 
back controllers that recover the aircraft from the angle-of-attack excursion, while 
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maintaining acceptable flying qualities, and without saturating the actuator, or ex- 
ceeding a certain actuator rate, empty? 

3. Proposed Numerical Solution 

In this section, we make the conjecture that the problem at hand is analogous 
to the example presented in the introductory section. The center-of-gravity location 
plays the role of the parameter influencing the stability of the open-loop plant. The 
tail volume plays the role of the parameter influencing how controllable the plant is 
from the input where feedback is to be applied. It is shown that when the plant 
parameters are fixed, the problem can be formulated in terms of an LMI, which is 
affine in both the controller parameters and actuator limits. Minimizing the actuator 
rate limit is a convex optimization problem that can be solved exactly. For one set of 
plant parameters, let the result of the process be characterized by the final values of 
all of the parameters. Repeat the process for different sets of values of the plant pa- 
rameters. Then, an assumption is made that the points lie on a smooth hypersurface. 
If the grid of plant parameters is sufficiently fine, the shape of the hypersurface can 
be discerned. The projection of the hypersurface into the three dimensional space 
spanned by tail volume, center-of-gravity location, and actuator rate limit defines a 
surface in the Tail Sizing Design Space that is extremely useful in the design of the 
aircraft. 

Let 



$1 {Ql {Zo, X cgo , Vh 0 , Vc 0 ) , Uh 

max ? Uh max ,u 

Cmax ? 'U'Cmax ? v o) = 

{W,Y > 0 : 

Y Y T [0C Tc 0] t + W T [D rc 0] T 

. [0 CV c 0] Y + [D rc 0] W u c 2 _ 

1 i% 

> 0 , 

vo Y 
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where 



Let 



Y 

[0 0 Cr h ] Y + [0 Dr J W 



Y T [0 0 C r f + W T [0 D r h f 



Ut 



> 0 , 



1 v% 
V 0 Y 



> 0 , 



(sin 0)(^V' + BW) + (sin <t>){AY + B\V) T -(cos <*>)( .4 y + BlV) + (cos <fi)(AY + M) T 
(cos^)(^r + BW) - (cos <f>)(AY + BlV) T (sin <t>)(AY + + (sin <t>)(AY + BW) T 

0 0 



+ m + (^y + Bvr) r + 2/3y 



< 0 , 



Y Y T [0 C ac Of 



[o c ac o ] y u 



Cmax 



Y Y T [0 0 C ah ] T 
[0 0 C ah ] Y u\ max 



> 0 , 



> 0. 



A = 



A B c C ac B h C\ 



a h 



0 A 



a c 



o 



0 0 



A 



a-h 



,5 = 



B a , 0 



0 B, 



a/i 



$ 2 {Qi (Zo, x cgo V}f 0 , V Co ), u h max ? ^hmax 7 ^ Cmax 7 ^ Cmax ? ^o) = 

jj4yt, Bk, Ck , /?, S' : 



* / cl cl 

ISO 

C a c C k o ul max 

R i cicl 

ISO 

CaA 0 u\ 



> 0 , 



> 0 , 



(11.73) 



(11.74) 
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where 



R I 

I S 

[0 C Tc A a }R + Cs c B a C k [0C rc ]A a 

R I 

I S 

[0 C re A a ]R + C Se B a C k [0 CjA a 



ClBlCl + R T [AlCj c 0] 

A T a [Cl 0 ] 



u 



Cmai 



CJBICI + R T [AlCl 0 ] 

A T ACT 0] 



HI 



l 



Xn 



xtS T 



xo R 

Sx 0 / 



I 

s 



> 0 , 



(sin<^)(n T + n) (cos<^)(-n + n T ) 
(cos <f>)(U - n T ) (sin<^)(II T + n) 



o 

o 



n T + n + 



R 

i 



i 

s 



2/3 



> 0 , 



> 0 , 



< 0 , 



(11.75) 







A B c Ca c B h Ca h 




R _ 0 






A 


BcCa c 


B h Ca h 










0 Aa c 0 


R + 


°CLc u 

0 Bq. l 






0 


Aa c 


0 






n — 




0 0 




L a h J 






0 


0 


Aa h 






1 1 — 














A 


B c Ca c 


B h C <*h 
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Then, a sufficient condition for the existence of a dynamic, output feedbac- 
controller, or static, state-feedback controller, that stabilizes the feedback system 
Qi ( Z 0 ,X cgo , Vh 0 , kc 0 )? does not exceed actuator amplitude and rate limits, and results 
in acceptable flying qualities in response to the angle-of-attack excursion defined by 
Vo, is for the sets, 4>j, or 4*2 to be non-empty. At this point, we observed that the 
same numerical results were obtained using a slightly different approach. Instead of 
minimzing a parameter u max in the LMI decision vector, we exploited the fact that 
the short period pole becomes more unstable as the center of gravity is moved aft. 
This allowed the use of a feasibilty algorithm vice minimization algorithm to acheive 



28 



the same results. The benefit was a significant, reduction in computational time. Now, 
the Plant Controller Optimization (PCO) problem considered in this section can be 
stated as follows: 



Maximize {x cg } 

Subject to: 

U 0 ,x cg , Vh 0 , Vc 0 ) = 0, 

H(Xo,U 0 ) = Z 0 , 

{Y,W) £ $\{Gl (x C g,VH 0 ,Vc 0 ),Uh. max ,Uh. max , ^ Cmax 7 ^Cmax v o), (H-77) 

or 

Ckt R) S) G ^2 {Ql Hoi^Co) ’fV'hmax '> U h m ax’> ^ Cmax 5 ^Cmax v o)- 

(11.78) 

A solution to this PCO problem includes a linear controller that stabilizes the feedback 
system Gi (Zo, x cgo , V// 0 , Yc 0 )-> and meets actuator limit requirements, as well as a 
providing a maximum aft center-of- gravity location. 

The numerical solution used to map the Tail Sizing Design Space involved a 
binary search over center-of-gravity stations: 

1. Fix V h ,V c ,Vt, 7- Let x cgmax = 1 and x cgm , n = 0. 

2. X cgo = (x c g max + X cgm,n)/ 2- 

3. Solve ^(Ao, Uq , x cgo ) = 0, 7"f(Ao, Uq ) = Zq for Ao, Uq. 

4. Obtain A(X 0 ,U 0 ,x cgo ), B(X 0 ,U 0 ,x cgo ) and form Gi (Z 0 , x cgo ). 

5. Solve for (K, W) € <M Gi (x cg ), u hmax , u hmax , u Cmax , ii Cmax , v 0 ) or, 

6. Solve for (A k , B k , C k , R, S) <E 4*2 (Gi (x cg ),u hmax ,u hmax , ^ Cmax 5 ^ Cmax "> Vo)- 

7. If no such (V, IT) exist 

• X eg max ~ x cg 0 -> 

else 
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• x cg m ,n — x cg 0 - 

8. If X C g max X Cgmxn ^ S° 2. 

9. Increment Vh, go to 1. 

4. Results - Rigid Body HSCT 

This design methodology was applied to a number of aerodynamic models 
representative of current high-speed civil transport designs. The results are shown 
for the aerodynamic model termed Ref A (see Appendix A). Appendix A details how 
a convenient build-up of the nonlinear longitudinal dynamics in terms of wing-body 
and tail contributions facilitates the process of mapping the Tail Sizing Design Space. 

The state feedback synthesis LMIs (set 4>i) were used. Figure 8 shows an 
optimization run for a single tail volume of 0.2, no canard, and peak actuator rate 
limit of 30 degrees per second. As the center-of-gravity location is moved aft, the peak 
actuator rate approaches the limit. The actuator amplitude remained well below its 
saturation limit, which was set at 20 degrees of travel from trim. 
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Figure 8. Fixed Tail Volume and Peak Actuator Rate Limit 

Figure 9 details the process for one tail volume and a sweep of peak actuator 
rates. Again, the tail volume was 0.2 and the peak actuator rate limit was incremented 
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from 5 to 30 degrees per second. Notice, initially the center of gravity moves quickly 
aft with only small increases in the peak actuator rate required. However, at some 
point, the peak actuator rate required becomes very sensitive to changes in the center- 
of-gravity station. For the Ref A model, this break point is in the vicinity of 0.565 
percent mean aerodynamic cord. 
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Figure 9. Fixed Tail Volume, Sweep of Peak Actuator Rates 



This process is repeated for the range of tail volumes of interest. Figure 10 
shows the results of this process for a tail volume of 0.1, and a tail volume of 0.2. 
A two dimensional representation of the data is cumbersome. When repeated for 
numerous tail volumes, the Tail Sizing Design Space could be viewed in three dimen- 
sions. Figure 11 shows the result of fitting a surface to data obtained for a range of 
tail volumes from 0.1 to 0.3, and a range of peak actuator rates from 5 to 30 degrees 
per second. The surface represents a lower bound on the peak actuator rate required 
by the feedback control system to recover from the angle-of- attack excursion, for var- 
ious combinations of center-of-gravity location and tail volume. An upper bound in 
the Tail Sizing Design Space would be a fixed limit on the peak actuator rate avail- 
able. Figure 12 shows the inclusion of this plane for a value of 20 degrees per second. 
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The volume below the plane and above the curved surface represents the Tail Sizing 
Design Space , where linear state-feedback controllers are known to exist that meet 
design requirements. Figure 13 shows the intersection of the plane representing the 
peak actuator rate available of 15 degrees per second, with the curved surface rep- 
resenting the minimum peak actuator rate required for a rigid-body HSCT with no 
canard. This is seen to be an aft center of gravity limit line on the standard center 
of gravity versus tail volume (scissors) plot. 
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Figure 10. Two Tail Volumes, Sweep of Peak Actuator Rates 

At this point, the canard volume was fixed at 0.05. Following the procedure 
outlined earlier, the lower surface of the Tail Sizing Design Space was mapped for 
a range of horizontal tail volumes from 0.1 to 0.3, and for a range of peak actuator 
rates from 5 to 40 degrees per second. Amplitude and actuator rate limits for the 
canard and horizontal tail were matched. Figure 14 shows the resulting lower bound 
in the design space. Figure 15 shows the results for Ref A HSCT with and without 
the canard. Figure 16 shows a slice of the two surfaces in Figure 15 at a peak 
actuator rate limit of 15 degrees per second. This should be familiar as a conventional 
scissors plot. Of course, as one would expect, the design space increases with the 
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tail volume 



Figure 11. 3D Tail Sizing Design Space 

additional control power available from the canard. The pertinent question is whether 
or not there is any benefit in spreading the control power for and aft, so to speak, 
or if the Tail Sizing Design Space would have increased just as much had the tail 
volume alone been increased by 0.05, and the canard not added. In general, it makes 




Figure 12. 3D Tail Sizing Design Space with Peak Actuator Rate Limit 
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Figure 13. 2D Slice of the Tail Sizing Design Space 

more sense to compare the competing configurations in terms of equal amounts of 
combined horizontal tail and canard surface area. For instance, profile drag is more 
closely related to the surface area of the control surfaces, among other things, and 
the design goal might be to minimize drag for the same aft center-of-gravity station 
and actuator rate limit. For this example, which utilized the Ref A data, the distance 
from the vehicle’s wing-body neutral point to the aerodynamic center of the canard 
or horizontal tail was the same. Therefore, comparisons in terms of normalized area 
or volume are equivalent. 

Figure 17 compares the two configurations, the first without a canard and the 
second with a canard. The percent change in total control volume required in going 
from a configuration without a canard to a configuration with a canard is shown as a 
function of the aft center-of-gravity limit. The same flying quality requirements and 
actuator limits were used. For the rigid-body HSCT model, the benefit gained from 
the inclusion of a canard is about a 10 percent savings in total control volume. 

This process was repeated, using the output feedback synthesis LMIs (set <J> 2 ). 
Figure 18 shows the resulting Tail Sizing Design Space obtained utilizing dynamic, 
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Figure 14. 3D Tail Sizing Design Space with Canard 



full-state, feedback controllers. Figure 19 compares these results with those obtained 
utilizing static, state-feedback controllers. The intersection of the two surfaces with 
the plane representing a target rate limit of 25 degrees per second is shown in Fig- 
ure 20. By inspection, it is clear that the limits of the Design Space are nearly 




Figure 15. Tail Sizing Design Space With and Without a Canard 
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Scissor Plot 




Figure 16. 2D Slice of Tail Sizing Design Space With and Without a Canard 

identical. This is an important and somewhat surprising result considering the added 
complexity of the output-feedback controller. 




Figure 17. Decrease In Total Control Volume Through the Addition of a Canard 
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Figure 18. 3D Tail Sizing Design Space - Dynamic Controller 



5. Problem Formulation - Aeroelastic Model 

The development of an integrated aeroelastic aerodynamic model is well docu- 
mented in work by Waszak and Schmidt [Ref. 45]. The nonlinear equations of motion 
were derived using a Lagrangian approach and some standard simplifying assump- 
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Figure 19. Tail Sizing Design Space Comparison: Static versus Dynamic Controllers 
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Figure 20. 2D Comparsion: Static versus Dynamic Controller 



tions. The elastic deformation is assumed sufficiently small such that linear elastic 
theory holds. Informally, the flexible deformation is captured by the contribution of 
an infinite number of states termed the generalized elastic coordinates. Associated 
with each generalized elastic coordinate is an in vacuo mode shape. In laymens terms, 
this describes the shape of a single elastic mode when fully deflected. Of course, for 
practical purposes, only a finite number of elastic modes are retained for analysis. In 
a body-fixed reference coordinate system, the resulting elastic airplane equations of 
motion for the longitudinal dynamics are: 



m(u + qw + gsinO ) 


— Qxi 


m(w — qu — gcosd ) 


= Qz, 


lyyQ 


= Qe , 


mi(q t + 2(u n rn + ufai) 


= Qi hi 



( 11 . 79 ) 



where, 




Q j = generalized forces , 

im = generalized mass , (11.80) 

77 ; = generalized elastic coordinates. 

Let X, and Z be the total aerodynamic and propulsive forces along each of the 
axis in the body-fixed reference frame. Similarily, let M be the total aerodynamic 
and propulsive moment about the body-fixed y-axis. Define 8x , 8z to be virtual 
displacements along the x and 2 axis, and define 80 be a virtual rotation about 
the y axis. When combined with the virtual pertubations of the elastic generalized 
coordinates, the virtual displacements are called the generalized coordinates. The 
Principle of Virtual Work is used to expand the generalized force terms in (11.80) 
through the use of the generalized coordinates. It states that 



Qi = 



d{8W) 



(11.81) 



d(S qi ) ’ 

where SW is the work associated with an arbitrary virtual displacement of the gen- 
eralized coordinates, Sq{. This virtual work done by the aerodynamic and propulsive 
forces, relative to the inertial reference frame, expressed in a coordinate system at- 
tached to the body of the aircraft, is: 



SW 



r ~ 00 

= XSx + ZSz + \M + ( zX — xZ ) 89 b + / P{x , z) • ^ <f) t 8rndS. (11.82) 

J s i= i 



The last term represents the work done by the distributed surface pressure 
due to virtual displacements of all of the elastic generalized coordinates. 

Applying (11.82) to the expression for virtual work, the generalized forces are 
seen to be 



, _ d(*W) _ 

x d{8x) 
c>(<5W) _ 
2 d(8z) 



(11.83) 

(11.84) 



39 



(11.85) 



Q - = Wri = a(|o ( / s P(I ’ Z) • (IL86) 

At this point, a method needs to be chosen to determine the aerodynamic and 
propulsive forces and moments. We assume that a knowledge of wing-body, tail, and 
canard aerodynamic stability and control derivatives for the rigid-body aircraft are 
available. Furthermore, the vehicle is assumed to have one flexible mode, with its 
mode shape and generalized modal mass known. Generally, the first symmetric mode 
shapes of all the HSCT designs are very similar. Each mode increases the number of 
states in the linear model by two, which subsequently increases computational times 
required by the Tail Sizing Design Tool. Since the modes retained in the linear model 
are those that we wish to actively control, the maximum number of modes retained 
would be three or four. This method can be combined with other methods, [Ref. 45], 
and used to capture the residuals from the truncated modes. The method used here 
works well for capturing the interaction of the rigid-body states with the flexible- 
body states for the first few symmetric modes. It is worth noting that the LMI based 
Tail Sizing Design Tool tool can utilize any method that is capable of generating a 
linear aeroelastic model at a specified center of gravity, tail volume/canard volume, 
and flight condition. The main contribution of this method is its suitability to the 
iterative nature of the numerical solution, and its use of widely available rigid-body 
aerodynamic stability and control derivatives. 

Therefore, the aerodynamic and propulsive forces are expressed in terms of a 
body-fixed reference frame as 

X = L sin a — D cos a + T x , 

Z = —Lcosa—Dsma + T Z i (11.87) 

where L and D are the lift and drag aerodynamic forces, a is the angle of attack, and 
Ti is the component of thrust in the i t k direction. The total lift of the airplane is 
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L — L c + L w b + L t , 



( 11 . 88 ) 



where the subscripts are meant to be suggestive of the canard, wing-body, and tail 
contributions. In coefficient form, the expression becomes, 

Cl = I^Clc + Ch wb + ~gC i t , (11.89) 

where 5, S tj and S c are the reference area of the wing, tail and canard respectively. 
For the problem at hand, (11.89) is rewritten explicitly in terms of a tail and canard 
volume. Let c denote the reference mean aerodynamic cord, l t denote the distance 
from the aerodynamic center of the tail to the aerodynamic center of the wing-body, 
and l c be a similar length associated with the canard. Then, (If. 89) can be expressed 
in terms of control volumes as 



Cl = V C jC Lc + C Lwb + V H jC Lt . (11.90) 

l c H 

The individual lift coefficients are expanded in a first order Taylor series ex- 
pansion about the components of the state vector and control inputs. Local changes 
in the free stream dynamic pressure have been incorporated into the local lift-curve 
slope derivatives. For example, the coefficient of lift for the tail is expressed as follows: 

C L t = Cl 1o + Cli^ y + Ci ta a + CLt^& + Ci tq q 

oo oo 

+ y. C LtVt i n + Cl,,- Vi + Ci tuh Uh. (11.91) 

i=i ' t=i 

Most of the terms in (11.91) are familiar as aerodynamic stability derivatives 
for a rigid-body aircraft. Some, such as Cx, , may be new to the reader. In order to 
derive an approximate expression for these aeroelastic stability derivatives, we need 
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to introduce the concept of a mode shape. It is assumed that the general elastic defor- 
mation of the unconstrained body can be described as the product of two functions, 
one a function of only spatial coordinates, and the other a time dependent function. 
Informally, the mode shapes, or free vibration modes, are the spatial function, and 
the generalized elastic coordinates are the time function. Then, the local relative elas- 
tic displacement is described in terms of the mode shape, and the generalized 

coordinate, 7?,(t), as 



OO 



d = ^Z<t>i{ x ) ■ »?*(*)• 

1=1 

Similarly, the local relative elastic torsion is described as 



(11.92) 



y- d4>i( x ) /j\ 



1 = 1 



(11.93) 



Expressions 11.92 and 11.93 can be used to obtain formulae for the force and 



moment dependence on flexible motion. Again, considering the tail contribution to 
lift as an example, the lift-curve slope dependence at the tail on ij is computed from 
the rigid-body aerodynamic stability and control derivatives, and from expressions 



11.92 and 11.93 as 



d<j>i{x) 

U L trJx — ^L ta ( ^ )x act 1 

r -r 

L ' Lt n, ~ ° L ‘a V y 



(11.94) 

(11.95) 



where x aCt is the normalized location of the mean aerodynamic center of the tail. The 
canard and wing-body terms are treated in an analogous fashion. Once the total lift 
is calculated, an assumed knowledge of a drag polar is used to calculate the total 
drag. 

Similar to the build-up of the total lift on the vehicle, the total pitching mo- 
ment on the vehicle is expressed in coefficient form. Again, a first order Taylor series 
expansion about perturbations of the state variables and control inputs results in 
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(11.96) 



Cm — Cm 0 + CM a a -f Cm a,® + Ca/,<7 + CM Uh uh 

oo oo 

+ Cm uc U c + Ca/„, Vi + 7 A 

1=1 1=1 

The individual derivatives are written in terms of the center-of-gravity location 
( x cg ), aerodynamic center of the wing-body (x aCu;(> ), and familiar control volume terms 
[Ref. 14]. For instance, 



Cmo, = C La (x cg - x aCwb ) - V H 

+ V C M±, + 



ac L , 

da 



da 



da ’ 



(11.97) 



and 



Cm„, = C L ,Mc-x. c , b )-V H dCL ‘ 



+Vc 9 -§± + 

OT)i drji 



dr]i 



(11.98) 



which are composed of stability and control derivatives that were either assumed 
provided or previously computed. 

The remaining term to consider involves expansion of the generalized force 
associated with the elastic generalized coordinates. Recall, 



Q„ = -^-r(J s P(x,z)-}iS V ,dS). ( 11 . 99 ) 

As a first step, the integral expression in (11.99) is separated into three parts: 



f P(x, z) • (j)iSrjidS = / P(x, z) ■ (friSrjidS + f P(x, z) ■ (f>iSriidS 

JS Jc Jwb 

+ jp(x,z)-4i6iudS. ( 11 . 100 ) 
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The pressure distribution over each surface is approximated by a point force acting at 
the aerodynamic center of the lifting surfaces of the aircraft, and the axial component 
is ignored. Therefore, the three integrals above are approximated as 



/ P{x,z) ■ faSriidS « Fz c M X aCc)t>Vi + Fz wb M X ac wb ){>Vi + F Z ,M X ac t )6m- 

JS 

(11.101) 

Then, using (11.99), a reasonable approximation of the generalized forces is 

Qvt ~ F > Z c 4 > i{ x aCc ) “b P 1 Z w b4 > i{ x ac wb ) “h Fz t <f>i( x ac t ) • (11.102) 

Each term in (11.102) is expanded in a first order Taylor series expansion. The force 
on the tail, for instance, becomes 



Fz t = Fz t ^ y + F Zta a + F Zta a + F Ztq q + F Ztg 6 

OO OO 

+ ^2 F Ztt)t 1 H + ^2 F Zt r/i + F Ztuh u h . 

i = 1 t=l ' 



(11.103) 



Recall, the generalized forces, Q v , drive the dynamics of the elastic generalized coor- 
dinates. The rigid-body states couple into the elastic states via terms like 



— - c — 

F Zta 4>iOt{x aCt ) = cos aqSVHjC Lta <}>i(x a c t )a, 

n 

and the control inputs couple into the elastic states through terms like 

— - c - 

F ZtUh <t>i{x act )u h = cosaqSVnjC LtUh <i>i{x a c t )u h . 

The remaining terms couple the elastic states among themselves, 
fifth term in (11.103) for an example, we obtain 



( 11 . 104 ) 



( 11 . 105 ) 



Using the 
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(11.106) 



F ZtVi = cos aqSV H jC Lta ( ) JaC( ^i( J ac), 



and 



= cos a q SV H yC Ll (11.107) 

T, » it a v 

This completes the description of the modeling of the nonlinear aeroelastic 
dynamics. The intent was not to document each term. That level of detail is left to 
Appendix A. Instead, the flavor of how to incorporate the influence of a few elastic 
modes was demonstrated. The method assumed a knowledge of wing-body, tail and 
canard aerodynamic stability and control derivatives for the rigid-body, and a set of 
mode shapes and generalized modal masses. 

The nonlinear aeroelastic equations of motion were linearized at an equilibrium 
point determined by the flight condition, canard volume (possibly zero), tail volume, 
and center of gravity. The resulting linear model is given by: 



where 
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(11.108) 



x R = 
x E = 



T T 



u w q 9 
iT 

V V 



Generalized elastic coordinates are somewhat lacking in physical intuition, and 
are not directly measured by any sensor suite. The relationships between sensed pitch 
angle and pitch rate at a given local body station (subscript /), the rigid-body states 
(subscript R ), and the generalized elastic coordinates are: 



9i 



J )xiVt, 

i = i ux 



(11.109) 
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( 11 . 110 ) 
(II. Ill) 



qi = Oi, 

~ 6r dx )xiVl ' 



Let 



T = 



BxA ®Ax2 

0 0 1 0 0 (§)„ 

0 0 0 1 (§)„ 0 



( 11 . 112 ) 



define a similarity transformation that replaces the generalized elastic coordinates, 
77 , 77 , with a more physically meaningful pair of states, such as qt and 6[. The location 
chosen for the placement of the local pitch rate and pitch angle sensors was the cockpit 
station. Typically, an area is sought where the mode shape slope has its most positive 
value [Ref. 1 ]. The linearized aeroelastic model now becomes 
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(11.113) 



Using Ref A data from Appendix A, the behavior of a rigid-body HSCT aircraft 
is compared with an aeroelastic model with the first elastic mode retained. Table I 
compares the eigenvalues of the two models. Figure 21 highlights the undesirable 
effect of the elastic motion. It shows the pitch rate sensed at the cockpit to a step 
input of the elevator for the two models. 

The frequency separation between the elastic mode and the short period dy- 
namics is approximately 1 Hertz, and the damping of the flexible mode is only 0.02. 
Typically, attempts are made to attenuate the feedback prior to excitation of the 
flexible dynamics. On large transport aircraft with the flexible dynamics close in 
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Short Period 


Long Period 


Flexible Mode 




frequency 


damping 


frequency 


damping 


frequency 


damping 


Aeroelastic 


0.90 


.75 


0.14 


0.55 


6.7 


0.02 


Rigid-body 


0.91 


.75 


0.14 


0.07 


n/a 


n/a 



Table I. Eigenvalues of Aeroelastic and Rigid-body Models 




Figure 21. Pitch Rate at Cockpit, Aeroelastic vs. Rigid Body 

frequency to the short period dynamics, this is hardly possible. Furthermore, even 
if suitable notch or low-pass filtering within the control loop could be attained, the 
extremely light damping of the flexible modes results in problems in terms of gust- 
induced structural responses and fatigue life. Therefore, the problem posed will be 
one of actively controlling the flexible modes retained. Generally, this will entail im- 
proving the damping of the flexible dynamics, while ensuring stability of the short 
period dynamics as the center of gravity is moved aft. 

As before, actuator dynamics are appended to the aeroelastic model. The 
Plant- Controller Optimization problem formulation is exactly the same as discussed 
in section 1. 
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6. Results - Aeroelastic Model 

The results are shown for the aeroelastic dynamic model termed Ref B in 
Appendix A. Ref B utilizes the same rigid-body stability and control derivatives as 
Ref A. The effect of a single, symmetric, flexible mode was incorporated. The mode 
shape is shown in Figure 22. It was experimentally determined that requiring the 
feedback controller to increase the damping of the flexible mode to 0.05 resulted 
in acceptable damping of the short period mode. Once again, the LMI based Tail 
Sizing Design Tool was used to map out a surface in the Tail Sizing Design Space. 
Figure 23 shows the design space for Ref B without a canard. An upper bound in 
the Tail Sizing Design Space is shown by the level plane at a peak actuator rate of 
25 degrees per second. Figure 24 compares the aeroelastic model to the rigid-body 
model. Obviously, considerably more actuator rate is required for the aeroelastic 
model in order to recover from the angle-of-attack excursion, while actively controlling 
the flexible dynamics. Intuitively, this seems obvious. The Tail Sizing Design Tool 
provides a metric for comparison. For example, consider a target tail volume of 
0.20 and center-of-gravity station of 45 percent mean aerodynamic cord. This is well 
within the Tail Sizing Design Space of the rigid-body model, but outside that of the 
aeroelastic model. The aeroelastic model requires either an extra 0.05 increase in tail 
volume, or 14 degrees per second increase in peak actuator rate, to move within its 
Tail Sizing Design Space. In some respects, this is the price to be paid for actively 
controlling a non-rigid vehicle. 

Next, the effect of the addition of a canard to the aeroelastic model is explored. 
As before, a canard volume of 0.05 was selected, and the design space was determined. 
The rest of the design requirements remained unchanged, and the results are shown 
in Figure 25. Figure 26 compares the design space of the aeroelastic model with 
and without the canard. Figure 27 shows the resulting aft line on the tail volume 
sizing plot when the two are cut at a peak actuator rate of 25 degrees per second. 
Of course, the added control volume due to the canard allows a further aft center-of- 
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Figure 22. I st Symmetric Mode Shape 



gravity location for a similar actuator rate limit. 

The relative effect of adding a canard is addressed by comparing the total 
control volume, (Vh + Vc), required for a given aft center-of-gravity limit and peak 
actuator rate limit. The results are shown in Figure 28 as a percent reduction in 




c.g. (%c) 



tail volume 



Figure 23. 3D Tail Sizing Design Space: Aeroelastic Model 
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Figure 24. Tail Sizing Design Space: Aeroelastic versus Rigid Body Model 



control volume required in going from the configuration without a canard to the 
configuration with a canard. Since the lengths of the tail and canard moment arms 
are the same, Figure 28 could also represent a savings in total control effector surface 
area. Experience has shown that the use of canards is desirable for flexible aircraft 
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tail volume 



Figure 25. Tail Sizing Design Space With Canard 
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Figure 26. Tail Sizing Design Space Comparison: Canard On and Off 



[Ref. 1]. This method provides a metric to quantify that benefit. In this example, 
the inclusion of a canard is 300 percent more effective when added to the aeroelastic 
model then when added to the rigid-body model. 

As before, the process was repeated using the output feedback LMls. The 




Figure 27. 2D Comparison: Canard On and Off 
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Figure 28. Decrease In Total Control Volume Through the Addition of a Canard 

results utilizing the aeroelastic model support the conjecture that dynamic controllers 
of the same order as the plant do no better than static, state-feedback controllers. 
A representative comparison is shown in Figure 29. There the aft center-of-gravity 
stations are shown for a range of tail volumes and fixed canard volume. The peak 
actuator rate limit was 40 degrees per second. 

D. GUST RECOVERY 

A second dynamic requirement imposed on transport aircraft is the ability 
to recover from a severe gust. The rationale is that the open loop HSCT should 
not be so unstable that an unrealistically fast actuator and flight control system are 
required. The gust recovery criterion is similar to that of the high angle-of-attack 
excursion in many respects. Both are at their most adverse when the vehicle is at a 
slow speed. Also, the aft center-of-gravity configuration is limiting. However, unlike 
the high angle-of-attack excursion, which is assumed to occur essentially instantly, 
gust recovery criterion involves the vehicle penetrating a particular wind-shear over 
time. 
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Figure 29. 2D Tail Sizing Design Space Comparison: Static versus Dynamic Con- 
troller on an Aeroelastic Model 

The wind-shear occurs in the vertical plane, begins at zero velocity, builds to 
some peak value, and then dies out again. Since the shear is defined to occur as 
the vehicle crosses into a fixed region in space, the rate at which it is penetrated 
depends on the flight speed of the vehicle. In some industry documents, the gust 
profile is defined to reach its peak within 30 feet of travel of the vehicle. However, 
there is ambiguity as to the rate and peak value of the gust. The bottom line is that 
relatively late in the design/testing process, a flight simulation will be conducted. 
The vehicle will penetrate some gust, based on empirical weather data gathered from 
airports throughout the country. Failure at this stage involves a major redesign. 
What is required is the ability to quantify the concepts used in the definition of the 
design requirement, such as so unstable , and unrealistically fast actuator . The Tail 
Sizing Design Tool provides a bench mark to assess competing configurations that is 
directly related to recovery from a gust penetration. 

1. Additional LMI Considerations 

Recall the state feedback synthesis LTI system description. 
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x = Ax + Bu 
z = Cx + Du = < 
u = K x 

and suppose that the desired block structure of the feedback gain, K, is 



x — [A-\- BK]x 
z = [C + DK]x, 



( 11 . 114 ) 



K = 



K x 0 



( 11 . 115 ) 



A sufficient condition to assure the proper block structure of the feedback 
gain is to specify an appropriate block structure to the unknown variables in the 
corresponding LMIs. Recall, for the state synthesis feedback problem, the feedback 
gain is recovered as K = ITT -1 . Therefore, specify 



Y = 



Vi 0 

0 y 2 



( 11 . 116 ) 



such that 
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W 2 0 
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( 11 . 117 ) 



2. Problem Formulation 



The gust recovery requirement describes the face of the gust profile as a One 
Minus Cosine disturbance. An acceptable alternative to the One Minus Cosine profile 
is to fit two exponential functions to the profile. This captures the important rapid 
change in airmass velocity as the vehicle penetrates the shear. Unlike the One Minus 
Cosine disturbance, the gust velocity returns to zero after reaching its peak. The 
rate of decay of the gust velocity after reaching its peak, however, can be made quite 
slow. Therefore, the impact on the vehicle’s dynamics is minimal. 
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Let t p be the time it takes the HSCT to penetrate the shear, and let Gp be 
the peak vertical velocity of the shear. Then, the following functions define the two 
wind-shear profiles. 

One Minus Cosine Profile: 

/.(*) = f^l-cos(^)), 

= G„ 

Double Decaying Exponential Profile: 

h{t) = -C p exp (-Alt) +G p exp (-A2t) , 0 < t, (11.119) 

where 

Ai = 3 t p , and A 2 = (11.120) 

Both wind-shear profiles are shown in Figure 30. 



0 <t<t p (11.118) 

t > t p . 




Figure 30. Two Wind shear Profiles 

Let the vehicle’s linearized longitudinal dynamics be given by the LT1 system 
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( 11 . 121 ) 



v 



v 



Au -f* Bu , 

A u A w A g A e 




A ac t 
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as before. Based on (11.119), the state space description of the wind-shear is 
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( 11 . 122 ) 



where v g is the vertical velocity of the shear. Now, a model of the vehicle’s dynamics 
in an airmass whose vertical velocity is v g (t) is 



V = 



’ A 


(A u sin 0 O + A w cos 9 0 ) {A u sin 0 O + A w cos 0 o ) 
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V x n 



0 — G p G p 



The velocity states in (11.123) are inertially referenced quantities. Consider a 
similarity transformation based on the change of variables, 



C = Tri, 



where 
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C = TA/T-'C + TB IU , 

=: y4£ + Bu , 

Co = (11-123) 



Then, the first two states of ( are perturbations in true airspeed and angle of 
attack, both airmass relative quantities. This makes the first four states of ( natural 
quantities to be considered for feedback. However, for all practical purposes, the gust 
states are not measured, and should not be considered for feedback. This implies that 
the feedback considered should be of the form: 



I< = 



AT 0 



Ci 

C2 



(11.124) 



where ( 2 corresponds to the gust states, and Ci corresponds to all of the remaining 
states. 



3. Proposed Numerical Solution 

Let 
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(11.125) 



Now, replace $1 with $3 in (11.77), and proceed with the algorithm outlined 
in section (3). 



4. Gust Recovery Results 

This design tool was being developed concurrently with a NASA led multi- 
corporate effort to define a baseline configuration for an HSCT vehicle. The design 
engineers desired the ability to quickly asses the influence of changes in key parameters 
of the wind-shear on the final design. This necessitated that computational time be 
kept short, on the order an hour or less, since the intention was to investigate multiple 
scenarios throughout the day. 

In that light, the gust recovery criterion was applied to the numerical model, 
Ref A, with a fixed tail volume of 0.11 and no canard. Still to be finalized were the 
design requirements concerning the definition of the wind-shear. The design team 
desired to know how T sensitive their designs were to changes in peak values of the 
wind-shear and/or time to penetrate the wind-shear. The Tail Sizing Design Tool 
was used to compute aft center-of-gravity stations for differing wind shear profiles. 

Figure 31 shows the sensitivity of the aft center of gravity location to changes 
in the peak value of the gust. The penetration distance is one cord length. The peak 
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gust velocities were 30, 45, and 60 feet per second. The data indicates that at around 
20 degrees per second peak actuator rate, every 50 percent increase in the peak gust 
velocity pushes the aft center-of-gravity station forward 10 percent. 
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Figure 31. Increasing the wind- shear peak velocity moves the center of gravity limit 
forward. 



On the other hand, figure 32 shows the sensitivity of the aft center of gravity 
location to changes in the penetration distance. The peak gust velocity was 45 feet 
per second. A penetration distance of 1, l|, and 2 cord lengths is shown. Choosing 
an actuator rate limit of 20 degrees per second again, every 50 percent increase in 
wind-shear penetration distance pushes the aft center-of-gravity station aft 5 percent. 

On a percentage basis, the aft center of gravity limit is about twice as sensitive 
to changes in the peak wind-shear velocity as it is to changes in the penetration 
distance. 



E. CONCLUSIONS 

The sizing of the horizontal control surface(s) for HSCT is a difficult problem. 
Traditionally, aircraft definition has taken place apart from feedback considerations. 
Controller design was a derivative process, with little contribution to the aircraft 
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Figure 32. Increasing the wind-shear penetration distance moves the center-of-gravity 
limit aft. 



definition. Of course, the aircraft were statically stable, and feedback control w r as 
viewed as a significant enhancement but not flight critical. This is an interesting 
example of an application where traditional methods simply lack the necessary tools 
to adequately define the aircraft configuration. Clearly, the aircraft definition will 
come first. However, it must be done with tools that provide a quantifiable knowledge 
of the impact on the demands of the feedback control system. 

The most intuitive metric to select, capturing the demands of the feedback 
control system, is the peak actuator rate required. The inclusion of this metric 
adds an extra dimension to the tail sizing problem. The two dimensional tail sizing 
scissors plot is inadequate, and a natural extension is the Tail Sizing Design Space. 
A conventional scissors plot can be recovered by viewing the intersection of the Tail 
Sizing Design Space with a level plane. However, the three dimensional space allows 
the designer to see how sensitive the aft boundary is to changes in actuator rate. 
Clearly, there is value in knowing when small changes in the center-of-gravity location 
result in large changes in the actuator rate required. 

While not solvable analytically, new efficient algorithms make the problem 
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tractable numerically. A new design tool is provided that provides the capability 
to quickly determine the Tail Sizing Design Space for a given aircraft configuration. 
This chapter demonstrates that many tail sizing problems can be formulated as LMls, 
and exploit new interior point algorithms to obtain exact numerical solutions. The 
iterative nature of the design process requires a tool that can generate solutions in 
a timely manner. Using the Tail Sizing Design Too/, the user can make adjustments 
to the aircraft definition, and quantifiably asses the impact on the demands of the 
feedback control system. 

In the first half of the chapter, two fundamental changes in aircraft definition 
were explored. Their influence on the Tail Sizing Design Space for a representative 
model of an IISCT was quantfied. The first was the use of canards in addition to a 
horizontal tail. The second was the effect of simple, flexible motion. This resulted in 
the testing the matrix of basic aircraft configurations shown in table E. Additionally, 
the effect of using a static, state-feedback controllers, or dynamic, full-state, output- 
feedback controllers, was investigated for each element in the matrix. 



Rigid Body-Tail Only 


Rigid Body-Tail with Canard 


Flexible Body-Tail Only 


Flexible Body-Tail with Canard 



Table II. Matrix of Aircraft Definitions Tested 



Numerical results suggest that canards provide a small benefit for a rigid body 
HSCT. Their use is more effective when used on a flexible body HSCT. The metric 
used to asses their effectiveness was the change in total horizontal control volume. 
In this example, an aeroelastic model realized a reduction in total horizontal control 
volume of approximately 30 percent through the use of a canard. The rigid-body 
model realized a reduction in total horizontal control volume of approximately 10 
percent. The use of a dynamic, full-state, output-feedback controller did not improve 
the results. 
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In the second part of the chapter, a different dynamic constraint from the high 
angle-of-attack excursion was addressed. The dynamic constraint was the penetration 
of a severe vertical wind-shear. The definition of the wind-shear profile was open. 
Therefore, the effect of changes in the wind-shear profile on the Tail Sizing Design 
Space for a rigid-body HSCT was explored. In this example, a boundary of the Tail 
Sizing Design Space was found to be twice as sensitive to changes in the peak level of 
the gust then to changes in the distance to penetrate the gust. 
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III. INTEGRATED GUIDANCE AND 

CONTROL 

A. INTRODUCTION 

In a great number of envisioned mission scenarios, Autonomous Vehicles will 
be required to follow inertial reference trajectories accurately in 3-D space. See, for 
example, [Ref. 36, 9, 11] and the references therein. Similar requirements emerge 
from the recent work at NASA on the descent trajectory synthesis for air traffic 
control [Ref. 43]. To achieve this goal, the following systems must be designed and 
implemented on board autonomous vehicles (AVs): i) navigation , to provide estimates 
of linear and angular positions and velocities of the vehicle, ii) guidance , to process 
navigation/inertial reference trajectory data and output set-points for the vehicle’s 
(body) velocity and attitude, and iii) control , to generate the actuator signals that 
are required to drive the actual velocity and attitude of the vehicle to the values 
commanded by the guidance scheme. 

The advent of GPS (Global Positioning System) has afforded AV systems en- 
gineers a powerful new means of obtaining accurate navigation data that is required 
for precise tracking of given inertial trajectories. However, traditional guidance and 
control schemes used to steer the vehicle along such trajectories may prove inade- 
quate in the case where frequent heading changes are required, or in the presence of 
shifting wind [Ref. 35]. Traditionally, such systems are designed separately, using 
well established design methods for control, and simple strategies such as line of sight 
(LOS) for guidance. See [Ref. 20] and [Ref. 29] for interesting applications to un- 
derwater and air vehicles, respectively. During the design phase, the control system 
is usually designed with sufficiently large bandwidth to track the commands that are 
expected from the guidance system. However, since the two systems are effectively 
coupled, stability and adequate performance of the combined system about nominal 
trajectories are not guaranteed [Ref. 35]. In practice, this problem can be resolved by 
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judicious choice of guidance law parameters (such as the ’’visibility distance” in LOS 
strategy), based on extensive computer simulations. Even when stability is obtained, 
however, the resulting strategy leads to finite trajectory tracking errors, the magni- 
tude of which depends on the type of trajectory to be tracked (radius of curvature, 
vehicle’s desired speed, etc.) [Ref. 35]. 

This chapter proposes a new methodology for the design of guidance and con- 
trol systems for AVs, whereby the two systems are designed simultaneously. This 
methodology has two main advantages over traditional ones: i) the resulting trajec- 
tory tracking system achieves zero steady state tracking error about any trimming 
trajectory, and ii) the design methodology explicitly addresses the problem of stability 
of the combined guidance and control systems. 

Earlier work based on this approach can be found in [Ref. 15, 42, 25]. In par- 
ticular, the authors introduce a methodology for the design of controllers for UAVs 
to track inertial trajectories that are given in space and time coordinates. The tra- 
jectories considered are equilibrium (also known as trimming) trajectories of AVs, 
which are helices parameterized by the vehicle’s linear speed, yaw rate and flight 
path angle. Furthermore, in [Ref. 15, 42, 25] it is shown the linearization of the vehi- 
cle error dynamics and kinematics about any trimming trajectory is time-invariant. 
Thus, the problem of designing integrated guidance and control systems for AVs to 
accurately track trimming trajectories can be solved by using tools that borrow from 
gain scheduling control theory, particularly those reported in [Ref. 27]. Within the 
framework of [Ref. 27], the vehicle’s linear speed, yaw rate, and flight path angle play 
the role of scheduling variables that interpolate the parameters of linear controllers 
designed for a finite number of representative trimming trajectories. The results in- 
troduced in [Ref. 27] on the D-implementation of gain scheduled controllers can then 
be used to obtain a combined guidance/control system such that the properties of the 
linear designs are recovered locally, about each trimming trajectory. An interesting 
and very important consequence of the ^-implementation is that it leads naturally to 
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a controller structure where the only exogenous commands required are the desired 
linear inertial position and the yaw rate, thus avoiding the need to feedforward the 
trimming conditions for the remaining state variables and control inputs. 

However, the technique presented in [Ref. 15, 42, 25] has a shortcoming which 
may be of concern for UAVs and AUVs when tracking trimming trajectories in the 
presence of changing air and water mass. Since the controllers described in those 
references achieve accurate tracking of trajectories defined in terms of space and 
time coordinates, the relative speed of the vehicle with respect to the air cannot be 
controlled externally, its value being computed internally as a function of the tracking 
error. In practice, this may lead to unacceptable performance in the presence of winds, 
since the change in the airspeed of the vehicle may result in stall or structural damage. 

Clearly, eliminating the time coordinate in the trajectory definition and using 
the vehicle attitude to null out trajectory errors while maintaining constant airspeed 
should resolve this problem. A similar approach has been introduced in a number of 
publications on robot control. Of particular interest is the work reported in [Ref. 40], 
where the subject of path following control for wheeled robots is addressed. See also 
[Ref. 35] for a detailed analysis of the stability of an autonomous underwater vehicle 
about nominal trajectories in the horizontal plane. 

In this chapter, these ideas are formalized within the basic framework for 3-D 
trajectory tracking controller system design developed in [Ref. 15, 42, 25]. Using 
the concepts outlined in [Ref. 40], the linear position of an AV is given in terms 
of its location with respect to the closest point on a desired trajectory, together 
with the arc length of an imaginary curve traced along that trajectory. Tracking of a 
trimming trajectory by the vehicle at a fixed speed is then converted into the problem 
of driving a generalized error vector - which implicitly includes the distance to the 
trajectory - to zero. Moreover, it is shown that the linearization of the generalized 
error dynamics about the corresponding trimming path is time-invariant. Using these 
results, the problem of trajectory tracking is posed and solved in the framework 
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of gain scheduled control theory, leading to a new technique for integrated design 
of guidance and control systems for AVs. The chapter summarizes the resulting 
design methodology, and illustrates its application to the design and implementation 
of a nonlinear trajectory tracking controller for the UAV Frog [Ref. 25]. Numerical 
simulations using a full set of nonlinear equations of motion of the vehicle show the 
effectiveness of the proposed techniques. 

The subject of trajectory tracking has also been addressed in the literature 
on control of nonholonomic vehicles. In [Ref. 44], the problem of tracking a nominal 
trajectory by a nonlinear system is considered. The key idea includes linearizing the 
nonlinear system along the trajectory, then using the resulting time varying lineariza- 
tion to obtain a time varying, state-feedback controller that locally exponentially sta- 
bilizes the system along the trajectory. The paper includes examples of applications 
of the proposed technique to a mobile robot and a front wheel drive car. The nominal 
trajectories considered in [Ref. 44] are not restricted to be trimming trajectories. 
However, all the examples presented consider the case of trimming trajectories only. 
Another approach is used in [Ref. 19], where a tracking problem for a surface ma- 
rine vessel is considered. Here the authors use feedback linearization with dynamic 
extension to obtain a controller to track trajectories that consist of lines and arcs of 
circles (a special case of trimming trajectories in the plane). 

The solution to the trajectory tracking problem proposed in this chapter dif- 
fers considerably from the ones introduced in [Ref. 44, 19]. Here, the key idea is to 
reduce the problem to the design of a tracking controller for a linear-time invariant 
plant utilizing a simple nonlinear transformation that inverts the vehicle kinematics. 
This poses no robustness concerns since kinematics are usually well known, partic- 
ularly in the case of air and underwater vehicles. It is important to point out that 
the application of the nonlinear transformation results in a nonlinear plant, whose 
linearization along trimming trajectories is time-invariant. Once in the linear setting, 
the designer is free to choose his favorite control synthesis technique to achieve the 
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desired closed-loop performance and robustness. The chapter provides a simple al- 
gorithm for implementing the linear controller on the nonlinear plant such that the 
properties of linear controller are preserved along each trajectory. This is in contrast 
to the approach in [Ref. 44], where the problem is reduced to that of designing an 
exponentially stabilizing state-feedback controller for a linear time-varying system. 
This leads to a controller design that is problem specific and does not address the 
issues of performance and robustness. On the other hand, in [Ref. 19] the authors 
point out that extending to air vehicles the feedback linearization technique use for 
trajectory tracking of the surface craft is difficult due to the unstable zero dynamics 
that are characteristic of aircraft. 

The chapter is organized as follows. Section B develops the rigid body dy- 
namics and introduces appropriate notation. Section C introduces the error dynam- 
ics necessary to solve the problem. Section D formulates and solves the problem of 
integrated guidance and control of UAVs for the class of trajectories introduced in 
Section C. Section E describes an application to the integrated guidance and control 
of the UAV Frog . Finally, Section F contains the main conclusions. 

B, RIGID BODY DYNAMICS 

This section begins with a review of the equations governing the motion of a 
rigid body, and specifically of an unmanned air vehicle. Notation familiar in robotics 
[Ref. 10], but not commonly used in the field of aircraft dynamics, is introduced. It 
is hoped that the familiar environment of the equations of motion of a rigid body will 
make the reader comfortable with the notation. It will be used extensively later in 
developing an integrated guidance and control algorithm. 

Let {/} denote an inertial reference frame. For the purposes of this devel- 
opment, the rotation of the earth and its associated Coriolis’ force can be ignored. 
Thus, a local tangent plane reference frame will be considered an inertial reference 
frame. Let {B} denote a body coordinate frame that is fixed with respect to the body 
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of the vehicle. Finally, let {IT} denote a stability coordinate frame aligned with the 
velocity vector of the vehicle. 

Free vectors can be expressed in whatever reference frame is most convenient. 
The velocity of the vehicle with respect to {/} is a free vector. When resolved in 
{£?}, it will be described by 



V 



u 



tT 



V w 



(HU) 



Free vectors can be transformed from one frame to another via the rotation matrix 
describing the relative orientation of the two frames. For example, it is common to 
describe the orientation of {£?} with respect to {/}, by the Euler angles ( <f > , 0 , ip). 
Let 



A 



<t> 

0 



then f R( A), or simply f R, is the rotation matrix from {/} to {B}, and given by 



B 

I 



R 



cos(0)cos(xp) cos(e) sin(V') — sin(0) 

sin(<£) sin(6?) cos( ip) — cos(<£) sin(^) sin( ip) sin(0) sin( V>) + cos(<£) cos(V') sin(<£) cos( theta) 
cos (<£) sin(e) cos( xp ) + sin(<£) sin( ip) cos(<£) sin(0) sin(pst) — sin(<£) sin(pst) cos(<£) cos(0) 



(II 1-2) 



As an example, consider the position P of the vehicle. Its velocity, ^P, in {/} can 
be expressed as 

d 

7t p = ° RV - 

Another common rotation matrix relates the orientation of {5} to the stability 
axis of the vehicle {IT}. The vehicle’s angle of attack, a, and side-slip angle, 0, 
define the orientation of {IT} with respect to {5}. The associated rotation matrix 
is termed ^ v R(a,0). As an example, if V is the velocity of the vehicle with respect 
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to {/}, resolved in {/?}, then 



0 

0 



= b’K(«,0)V 



is the same vector resolved in {VU}. 

The angular velocity of {B}, with respect to {/}, resolved in { B } will be 
represented by Q. It is the set of familiar body-fixed rotation rates, roll rate (p), 
pitch rate ( q ), and yaw rate (r). The body-fixed rotation rates are related to the 
Euler angles by 
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(III.3) 



Inverting (III. 3), the time rate of change of A is related to fl by the matrix Q(A), or 
simply Q, as 

id = ^ 

S(0) = /. (111.4) 



Using this notation, an application of Newton’s Law to the linear and angular 
motion of the vehicle yields 



T V 

dt 

T n 

dt 



F 

-n x u + — , 

m 

Jg'fi x J gH + ?b 1 N, 



(HI-5) 



where F and N are the total external force and moment acting on the vehicle resolved 
in {B}, m is the body mass, and 2g is the inertia tensor of the body resolved in {B}. 

Standard practice is to separate out the contribution of gravity, ( G ), from the 
force and moment terms. The remaining terms are expanded in a first order Taylor 



69 



series expansion about some choice of state variables and control inputs. A convenient 
choice of state variables is (u, /?, a, p, g, r), and standard control inputs are elevator 
(£e), aileron ( Sa ), rudder (Sr) and throttle (St). Applying this to (III. 5) results in 



±v 

dt 



-ttxV+f RG + Z-p\\V\\ 2 sK w {aJ) 
Zm 



+ 





Cd 0 




Cd u 


0 


CD a 




< 


Cy 0 


+ 


0 




0 






. ° L ° . 




Cl u 


0 


Cl q 





r 
















8t 




0 


c Dq 0 


ft + 


C D 6t 


Cd 6€ 


0 


0 




8e 




°Yp 


0 C\ r 


0 


0 


Cv 6a 


Cv 6r 




8a 


> 




0 


c Lq 0 




Cl ( , 


Cl 6 ' 


0 


0 
























8r 


J 



u 

fi 

a 



(hi. 



T n 

dt 





/ 

\p\\ v f 


sb 


0 


0 




/ 


-,|c 

o 

1 


x Xgfi -t- Xg* < 


0 


sc 


0 


7l w (at,/3) < 




Cm o 






0 


0 


sb 






Cn, o 



+ 



0 

C Mh 

0 



0 

Cm 6s 0 
0 Cn, 



Cu a C, 



6r 

0 

Cn 6t 





8t 


X 


> 




8e 




> 




8a 






- 


8r 









0 c h 0 




u 




Ci v 0 c lr 




6 0 0 

2||V|| U U 


+ 


Cm u 0 Cmc 




P 


+ 


0 c Mq 0 




0 s 0 

u 2||V|| u 




0 c h 0 




a 




1 

o 

i? 




1 

.ol; 

CN 

o 

o 

1 



n 



(III.7) 



70 



where 



s 

c 

m 

1b 

P 

b 

D, Y, L 
/, M, N 

C X0 

Cx, 



wing reference area, 
wing mean chord, 
mass, 

aircraft’s inertia tensor resolved in {£}, 
air density, 
wing span, 

normalized drag, lateral force and lift, 
normalized roll pitch and yaw moments, 
normalized nominal force or moment coefficient, 
stability derivative:^. 



Notice that some symbols in (III. 6) and (III. 7) have multiple meanings, and attention 
must be paid to the context in which they are used. The terminology is standard in 
aircraft dynamics. 

A convenient grouping of terms in (III. 6) and (III. 7) is to let 7i be a vector 
valued function of the control inputs, X be a vector valued function of all the terms 
affine in 7 Y, and T be a vector valued function of all the remaining terms. Then, 
equations III. 6 and III. 7 can be compactly expressed in terms of these functions. 
Subscripts are used to indicate the correspondence of the function to the appropriate 
linear or angular motion equation. As an example, 



Xv : = 



2m 






C D H 


C D Se 


0 


0 




0 


0 


C Y Sa 


cv. 


(III 8) 


Cl 6( 


Cl. 


0 


0 





Then, the complete equations of motion of the vehicle are expressed in this notation 
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as 



0 



±v =;F v (v,n,A) + iv(v,n)H(v,n,u), 
f t n =^h(v,ft, A ) + i Q (v,n)n(v,n,u), 

Tt P =bRV, 

It a = <2fl- 



C. GENERALIZED ERROR DYNAMICS 
1. Trimming Trajectories 



(HI 9) 



The objective of the guidance and control systems is to steer an autonomous 
vehicle along prescribed inertial trajectories Pc € P ? . Furthermore, we require that 
the vehicle be trimmed along any such trajectory. Let { C } define the coordinate 
system attached to the vehicle, and let A c define the desired inertial orientation of 
{ C }. The coordinate system { C } represents the desired inertial orientation of the 
vehicle along Pc- Therefore, at trim {B} = {C}. Next, we define the set £ of the 
trimming trajectories about which the vehicle is expected to operate: 



i Pc =h PVc , 

Ac = Q(Ac) ftc =: Qc^tc, 

Fv(Vc, &c, Ac) + lv(Vc, ClcffliVc, £lc, U c ) = 0, 
MVc,Slc, Ac) + Ia(Vc,ficMVc,n C ,U c ) = 0, 

(III.10) 

where Vc,ftc and U c denote the trimming values of V,n and U, respectively, and A c 
denotes the vector of Euler angles that describe the orientation of { C } with respect 
to {/}. 

From the definition of £ and equations III. 9, it can be concluded [Ref. 13] 
that trimming trajectories Pc £ £ correspond to helices: 




d 

dt 



A c 



0 

0 , 
Ipc 



(Ill'll) 
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(III.12) 



-v c cos( 7 c ) sin (ip c t) 
v c cos ( 7 C ) cos(ip c t) , 
-v c sin( 7 c ) 

where 

Tp c = desired turn rate, 
v c = desired inertial velocity, 

7 C = desired flight path angle. 




Thus, the trimming trajectories can be parameterized by the vector 

7/c = [n c ,^ c ,7c] T <5 T? 3 . (III. 13) 



In fact, given r/ c , we can determine the trimming values for Vc, ftc and Uc which will 
be important later for the controller synthesis. 

Integrating (III. 12) we obtain 



P c (()= in (4() 

— V c sin(7c)^ 

Equation III. 14 indicates that the radius of the helix is 

Vc cos(7c) 



be := 



V>c 



and the climb rate is 



(III. 14) 



(III. 15) 



K 



—v c sin( 7 c ). 



(III. 16) 



Using elementary ideas from differential geometry [Ref. 3], equation III. 14 can be 
reparameterized by considering the arclength s defined as 



s 



/ Wj t Pc(tWU 

/ v c dt. 



(III. 17) 
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Although the symbol for arclength is the same as that used in other sections for the 
Laplace transform variable, it is hoped that the context will make the distinction 
clear. Then, using (III. 15), (III. 16) and (III. 17), equation III. 14 is written as 



Pc{s) = 



b c cos(xp c ^) 

b c sm(7f c ) • 



(III. 18) 



Two parameters useful in describing the helix, Pc(s), are its curvature and its 
torsion. The curvature is defined in terms of the parameters in T] c as 

ac(s) = 

and the torsion as 

t(s) = 



2. Trimming Trajectory Coordinate System 

Now, let {v4} denote a Frenet frame attached to Pc(< s) [Ref. 3]. The x axis of 
the Frenet frame is termed T . It is a unit vector tangent to Pc{s) at s, and it points 
in the direction of motion on Pc{ s )- Its components are defined as 



i> c cos(^ c ) 



h c ij>, 



C 

2”’ 



sin 7 C 



(III. 19) 



(III. 20) 



T(s) 



dP(s) 

7s ’ 

-“*sin( 04 ) 

^ «*(&£) ’ 
he 

y c 

- cos (7c) sin(0c^) 

cos(7c)cos(V> c ^:) 

-sin(7c) 



(III. 21) 
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The y axis is termed N and is perpendicular to T. It is a unit vector pointing toward 
the center of the helix, and it is defined as 

dT(s) 

N(s) = ~dr /K{s) ’ 

-C° s (0 c £) 

-sin (V> c ^) • (III. 22) 

0 

The 2 axis is termed B , and is orthonormal to T and N. It is defined according to 
the right hand rule, B(s) = T(s) x N(s): 

sin(7 c )sin(V> c ^) 

5(5) = -sin(7 c )cos(V> c ^) • (Hl-23) 

cos(7 c ) 

The rotation matrix from {/l} to {/} can be defined by the vectors, T, N, and B as 

* A n = [T N £]. (III. 24) 

Let A Q be the angular velocity of {>1} with respect to {/}, resolved in {>!}. It can 

be shown that A Q, is equal to [rs 0 ks] t . Using A Q, and the Serret-Frenet Theorem 
[Ref. 3], the time rate of change of A 7Z is 

d(' A n) i(' A n ) : 



dt 



ds 



= [T N B] 



0 — k 0 

K 0 — T 

0 r 0 



= ' A KS( A n), 

where *S(T) is a skew symmetric matrix defined by [Ref. 10], 



(III. 25) 



S(Y) = Lx, 

0 — KS 0 

= ks 0 — ri 

0 ri 0 
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Let P(s 0 ) denote a point on the trajectory Pc € S. Define the error vector 



M e := fK(P - P(s 0 )). 



(III. 26) 



We will call Pc(sq) a projection of P onto Pc € £ when Me has the following form: 



M e = [ 0 y z] T . 



From the definition of Me, Pc(so) is also the point on Pc nearest to P. (For the 
discussion of when such a point can be uniquely determined, see [Ref. 40]). Let 
Pc(so) be the projection of P onto Pc € £■ Then, using the definitions of Me and 
{/!}, we obtain 



Therefore, 



7t p ~ » RV 

- it p ° + iA nM ^ 





s 




0 




0 


+' ns{ A M)M E + ! a n 


y 




0 




z 



(III.27) 



a b RV = fR' B RV, 





s 




0 




0 


+f r^pls^^Me +f R r A n 


y 




0 




z 
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(111.28) 



1 

o 

o 

1 

I 




s 


ZT 1 0 




y 


yr 0 1 




z 



a b 11V. (III. 29) 



s 




1 - yti 


0 


0 


y 


= 


ZT 


1 


0 


z 




yr 


0 


1 



3. Generalized Error Vector 

The design of an integrated guidance and control system for the plant Q in 
(III. 9) and the set £ involves obtaining linear models for Q along the trajectories in 
£. These models will necessarily be time- varying in the state space coordinates used 
in (III. 9). It turns out, however, that an appropriate coordinate system exists where 
the linearization of the plant Q along any trajectory Pc € £ is time-invariant. This 
coordinate system is discussed next. Let Pc, he € £ be given. Define 



V E = V — V c , 
11 E — il £i/C ) 
Pe = 



1 T 



(III. 30) 



s y z 

A e = Q~ l [A — Ac] , 

which can be interpreted as the generalized error vector between the vehicle state 
and the trajectory in S. Simple physical considerations show that the problem of 
following a trimming trajectory Pc (E £ at a fixed speed Vc is equivalent to driving 
the generalized error vector to zero. 

In order to have a complete description of the error dynamics, we need to derive 
an expression for the time rate of change of the generalized error vector, (III. 30). 
Expressions (III. 9) and (III. 29) suffice to describe the derivatives of the first three 
terms in (III. 30), since Vc and Qc are constant along trimming trajectories. The 
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fourth term, A#, is the only term that requires an explanation. Recall, 



Therefore, 



Since 



and 



we obtain 



A e = Q 1 [A — Ac] . 



37 a e = a c) + 37 C 1 [ a - Ac]- 



dt 



dt dt 



d 



dt 



s A - 



—Ac = Qc^tc, 



(III. 31) 



(III. 32) 



—A^ = f l — Q 1 Qcflc + -jjQ 1 [A — Ac] • 



(III. 33) 

Therefore, the generalized error dynamics are described by the system of equations 



Ge = 



j- t \ r E = ^Fv e {Ve^e-, Ab) + 1v e {Ve, GIe)'H{Ve- i flc, U), 

Tt = F^eWe-, ^E, Ae) + ^q e (Ve,He)'H(Ve,GIe,U), 



dp 

dV E — 



-I -1 



1 — 1/ AC 0 0 

ZT 10 

yr 0 1 



B 



nv, 



tAe = GIe — GIc — Q 1 Qcfic + C 1 QA.e, 



where 



(III. 34) 



Fv e {Ve, GIe, A e) — Fv(Ve + Vci GIe + flc 5 Q^e + Ac), 

and similarly for %v B , J~q e and Zn E - 

An important consequence of this choice of error dynamics is that the lin- 
earization of equations III. 34 along any trajectory Pc & £ is time-invariant. In order 
to derive the linearization of (III. 34), we require the following identity. 
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Identity 1: Let B R be the rotation matrix from {/?} to {I}, and let Q satisfy ^-(A) = 
Qft. Then for any vector X in 7v 3 



Jf(sttV) = -' b RS(X)Q-\ 

and 

^(fRX) = S(?RX)Q~K 

Proof: First vve observe that [Ref. 10]: 

f t {' B R) =' B RS((l). 

Next, from the definition of a cross-product we get 

S(X)Y = -S{Y)X 

for any X, Y £ IZ 3 . Now, using equations III. 37, III. 38, and the fact that 
constant vector, we obtain 

j t CsRX) = 

= ' b rs(si)x, 

= - ! B RS(x)n. 

Next, observe that 

jfi RX ) = r^ RX ^ 

= ±{‘.*x)ea. 

Equation III. 35 now follows by comparing equations III. 39 and III. 40. 

To obtain equation 1 11.36 note, 

(: ?R r 1 = CbR), 



(III. 35) 



(III. 36) 



(III. 37) 



(III. 38) 
X is a 



(III. 39) 



(III. 40) 
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thus 






-? RjfiR)? R, 



= -s(n)fR. 



(III. 41) 



Therefore, 



It (fRX) = - s(n)fRX , 



= S{fRX)tl. 



(III. 42) 



Moreover, using the chain rule, it follows that 

= ± { fRX)± A, 

Equation III. 36 follows readily from equations III. 42 and III.43. 
Finally, let A = 0, then 



(III. 43) 



I B n=f n = q = i. 



Now, using equations III. 35 and III.36, we get 



= ^(?/W)| A= o, 
= S(fRX)Q~\ 
= S(X). 



(III. 44) 



Since any trajectory Pc £ £ defines a trim condition for equations III. 6 and 
III. 7, the linearization of (III. 6) and (HI. 7) is naturally time invariant. Introducing 
the following notation, 



•4? = ^[^(.)+Ix(. )«(.)], 

B z = JL[1 x (.)H(.)\, 
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the linearization of (III. 6) and (III. 7) is written as 



-8V e = Al6V E + J&6Sl E + A v K 8\ E + Bv6U, 

at 

—fiOr = , ^cn_ , AQ J 

dt 



;Sn E = A% 8 V E + A% 8 n E + A*8A E + B Q 8U. 



(III. 45) 



We need only to derive the expressions for ^<5 >Pe, and j^8 A E to show that they are 
time invariant. 

Let 

t -i 



M = 



1 - yn 
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0 






ZT 
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0 


> 




yr 
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1 






(1- 


!/«)' 


-1 
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0 


—zt( 1 


-y«) 1 l 


0 


-yr(l 


1 

o 


1 



(III.46) 



Then, from (III. 34), 



dt 



P E = M{P e ) a b K{ A)V 

= Minv. 



(III. 47) 



Furthermore 
>• ■ 



jL[Minv\\jp E + dL\M A B nv\\M e + ^\M A nv\\,sv. 



(III. 48) 



Expanding the first term in (III.48), we obtain 



d 

dP E 



[M e KV]\ c = 



— I A vv 

dpJ cc7lVc ' 



dM 



dP El \c 



0 0 0 
0 0 0 
0 0 0 



dM | 


dM | 


9Pe 2 ' c 


9Pb 3 |c J 



cKVc, 



0 0 
1. 0 0 
| e 0 0 



0 

3 - 2 r(l-y «) _1 | 
dz I 

0 



0 0 

C o 0 

0 0 



c^Vc, 
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0 


0 


0 


0 


0 


0 


I 

o 


0 


0 



= (oKVc), 



K$nv c )x o o 
o oo 
-T($nv c ) x 0 0 

0 K 0 

0 0 — T 

0 -r 0 



0 0 0 
-t(£kv c ) x 0 0 



0 



0 0 



(III. 49) 



where the subscript c implies that the preceding gradient is being evaluated along the 
trimming trajectory. The subscript x implies the ^-component of the vector, and 

- iT 



Pe = 



Pr , Pe 2 Pe 2 



Now, from the definition of A e, it follows that 

dA 



dA 



= Q- 



E 



Next, using equation III. 50 and Identity 1, we get 

ik [Mm = 



d 



and 



= (g^TOW])|e, 

= ^*’0'" 
= -( fTl' B KS(V)Q-'Q )\ c , 

= -ins(v c ), 



= in. 



By combining equations (III. 49 - III. 52), we obtain 



d r 

7t sp ° 



0 K 0 

= &nv c ) x 0 0 -r 

0 — r 0 

- %RS{V c )8Ae+%K 6V e . 



SPe 



(III. 50) 



(III. 51) 



(III. 52) 



(111. 53) 

(111. 54) 
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The final term, -^8\ E , is obtained by observing that, along the trajectories 
Pc £ £, TtQ~ X = 0- Using simple algebra, vve obtain 

7t SA ° = 7t^§b ME - ( " L55) 

Now, since -^Ac = [0 0 tpc } 1 and from the definition of Q~ x (III. 3) and (III. 2), 
it follows that 

Q-'^c = Ifc. 

= ?RjA c . (111.56) 



Finally, using Identity 1 and equations III. 50 and III. 56, we obtain: 



dA 



(Q-'Uc) 



dA 



dt dAi 



“ P n 7t A ^ 



= S(flZ—A c )Q~ 1 Q\c = S( ft c ). (III. 57) 



The desired expression for jiSAe now follows from expressions III. 55 and III. 57 as 

jM e = Sn E - S{Slc)6A E . (III. 58) 

dt , 

Summarizing, the linearization of equations III. 34 along any trajectory Pc £ £ 
is time-invariant and given by 



' i*v E 

P^E 



Qi ( Vc ) = < 



P? 



Ay8V E + A^8VL E + A v k 8A e + B V 8U, 
Ay8V E + A&n E + A u h 8A E + BtfU, 



(£KV C ) X 



0 K 0 
0 0 — T 

0 -r 0 



8Pe, 



£KS(Vc) 8A e +£ K8V e , 



(III. 59) 



j- t 8A E = 8£l E — S(ttc)8A E , 

where (cP-Vc)x is the ^-component of the vector qTZVc- In the next section, the 
symbol Qi will denote the set of all linear plants Qi ( r / c ) associated with the set £. 
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D. TRAJECTORY TRACKING CONTROL SYSTEM DE- 
SIGN AND IMPLEMENTATION 
1. Linear Controller Design 

In the previous section, we have shown that the linearization of the nonlin- 
ear system Qe about any trajectory in £ results in a time-invariant plant £?/ (r? c )- 
Therefore, associated with the set £ there is a family of linear plants, Qi . These can 
now be used to synthesize a tracking (possibly gain-scheduled) controller C , which is 
designed to operate over all the trajectories in £. 

A common approach to the development of such a controller C requires de- 
signing a family of linear controllers for a finite number of linear plants in Qi , and 
then interpolating between these controllers to achieve adequate performance for all 
the linearized plants in Qi . During real time operation, the controller parameters are 
updated as functions of a gain-scheduling variable q = h(V, 0, A, P, U, rj c ), where h is a 
C 1 function. For example, a typical gain-scheduling variable, for the case of aircraft, 
is dynamic pressure. In that case, q = |p||Vj| 2 , where p represents air density and is 
itself a function of aircraft height. 

In what follows, we restrict ourselves to the idealized case where the description 
of each controller for each plant £?/ ( r / c ) is available [Ref. 39]. Therefore, we assume 
that the first design step produces the set 

Ci {Ci (<7c),<7c = h(Vc,£lc, he, Pc,U c ,T] c )}, 

given by (see Figure 33): 



8E 



Ci (q c ) = < 




c2 



su 



[6u Sy 8z] t - [8v c 8 y c 8 z c ) T , 
A cl (q c ) 8 X ci + B c ,(q c )[ 8 V T 8 ft T 8 A t e ] t 
+B C 2 {q c ) 8 X C 2 + B c z(q c ) 8 E, 

8E , 

C cl (q c ) 8 X cl + V cl (q c )[ 8 V T 8 £l T 8 AIY 

J rP ) c2{ ( lc)hX C 2 + 'D c z{ ( lc)£E, 



( 111 . 60 ) 



84 



where 8X C \ € R Uc , 8X C 2 € R m and m = dim(/ 7 ). The vector 8X C 2 represents the 
integrator states of the controller C[ (q c ). The vector 8X C 1 represents the remaining 
states of the controller Ci ( q c )• The velocity and positon error is 



8 v 


T 


1 


0 


0 




8V 




0 


0 


0 


8y 


= 


0 


0 


0 




8tt 
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0 


1 


0 


8z 




0 


0 


0 




1 

On 

> 

ft 

1 




0 


0 


1 



=: C x 



8V 

s n 

8Ae 



+ C28PE1 



8 Pe , 



(III. 61 ) 



where 8v c , 8y c and 8z c are introduced to determine how fast the error states, 8v , 8y , 
and 8z, go to zero. We further assume that the parameters of the controller are C 1 
functions of q c . 




)C(q c ) = 



A c \ (?c) 


B C 1 (?c) 


^c2(?c) 


Bcz(qc) 


Cel (?c) 


Del (?c) 


D c 2 (?c) 


D&i^Qc) 



Figure 33 . Set of Linear Controllers 



The structure of the controller Ci ( q c ) has the following important feature. 
Suppose the closed-loop system consisting of Qi ( q c ) and C\ ( q c ) given by equations 
III. 59 and III. 60 is asymptotically stable. Then, for a given q c , the controller Ci ( q c ) 
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will ensure zero steady-state error to a step input for the variables in SE. This includes 
errors in the vehicle’s inertial velocity v, and in the deviations from Pc, z and y. Zero 
steady state errors are achieved by integrating 6E. This structure is typical of tracking 
controllers, since they are designed to drive errors between step changes in reference 
commands and the corresponding plant outputs to zero in steady state. Notice that 
the block IC(q c ) (see Figure 33) may itself contain additional integrators. 

2. Gain Scheduled Controller Design 

Next, the family of linear controllers Ci ( q c ) must be implemented on the non- 
linear plant Q defined in Section B. This problem has been addressed in [Ref. 27] 
for the general class of nonlinear plants and for tracking controllers with the same 
structure as C\ ( q c ). In [Ref. 27], the authors formulated a so-called controller imple- 
mentation problem which will be repeated here for the problem at hand. 

Let T (Qi ( q c ), Ci (<7c)) be the closed-loop linear system that results from con- 
necting Ci ( q c ) to Qi ( q c ), and denote by T(Qi ( q c ), Ci ( q c )) the corresponding matrix 
transfer function. Let T (Q , C ){q c ) be the nonlinear closed-loop system that consists 
of C and Q , and let , C )(^ c ) denote its linearization about Pc E £, and denote 
by Ti(Q , C )(q c ) the corresponding matrix transfer function. With this notation, the 
controller implementation problem applied to the integrated guidance and control 
problem considered in this chapter can be stated as follows: 

Controller implementation problem: “Find a gain scheduled controller 

C such that for each trajectory Pc € £ 

1. the feedback systems %(Q , C )(q c ) and T(Qi ( q c ), Ci ( q c )) have the same 
closed-loop eigenvalues. 

2. The closed-loop transfer functions T;(£/ , C )(<7 C ) and T{Qi ( q c ), Ci (q c )) are 
’’equal.” 

Given the set C/ of linear controllers for the set Qi of linearized plant models, 
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we propose the following structure for the gain scheduled controller C (see Figure 34): 

' [0 y z] T = f 7 Z(P - Pc(so)), 

E = [v - v c y - y c z - z c ] T , 
iX cl = A cl (q)X ci + B cl (q)[j- t V r j- t n T ({l-Q-'A c f]T 
C:=l +B c2 (q)E + B c3 (q)f t E , (III. 62 ) 

f t x c2 = C cl (q)X c ,+V cl (q)[f t V T %Q. T ( 1 n-Q-'A c f) T 
+'B>c 2 {q)E + V c3 (q)j- ( E, 

u = x c2 . 



Recall, P c (so) is the projection of P onto the helix Pc G S. Comparison of Figure 33 
and Figure 34 indicates that the structure of the gain scheduled controller is easily 
obtained from that of the linear controllers. 




Figure 34. Gain Scheduled Controller 



We now make the following assumptions: 

Al. Dim(A' c2 ) = dim (U) = dim (E) 

A 2. The matrix 

sI-Ad{q) B c2 (q) 
—C c \(q) V c2 (q) 
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has full rank at s = 0 for each Pc 6 8. 

A3. The matrix pair (A c i,C c i) is observable. 

Assumption A\ implies that the number of integrators is equal to the number of 
control inputs. This is necessary if the controller is to provide independent control of 
the errors E using the control inputs U. Assumption A 2 implies that the realization 
(A c ut3 c 2,C c i,T) C 2) has no transmission zeroes at the origin. Finally, assumption A3 
guarantees that the state X ci is zero along the trajectories in 8. 

The main result of this section is stated next. 

Theorem D.l Suppose assumptions A 1, A 2 hold. Then the gain scheduled controller 
C given by equations III. 62 solves the controller implementation problem, i.e., for each 
Pc € 8 the following properties hold: 

1. The feedback systems Ti(Q , C )(q c ) and T(Qi ( q c ), Ci {q c )) have the same 
closed-loop eigenvalues. 

2. The closed-loop matrix transfer functions TfQ , C )(<? c )( 5 ) and T(Qi ( q c ), Ci (<? c ))(s) 
are equal. 



Proof: In the proof, we set the controller matrices X> cl , to zero. This does not 
change the results but considerably simplifies the algebra. 

Let Pc G 8 be given. Let 



A, 
A 2 
A3 



< 




A\ 


Ay 




A% 
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-s(n c ) 




0 - 


4ns(v c ) 






0 K 0 

0 0 — r 

0 — r 0 
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13 

(£KV C ) X 

8 T 

Then, (III. 59) is expressed as 



B v 

Bn , 

0 

^-component of q'TZ Vc 

SV E 

SQe i 
SA e 



Gi (vc) 



\ &sp 



The error signal is expressed as 



ArfT + BSU, 
A 2 8T + A 3 8P. 



8E = C 1 8T + C 2 8P, 



which, when substituted into the linear controller (III. 60), results in 



(III. 63) 



(III. 64) 



Ci (q c ) = 



8X cl = A cl (q c )8X cl + B c 2 8X c2 + [B el {g c ) + B^C^ST 

+B c3 (q c )C 2 8P, 

&8X C 2 = C,8T+C 2 8P, 



(III. 65) 



8U — C c \ {q c )8X c \ V c2 {q c )8X c2 . 

Consider the feedback interconnection of the linear plant (III. 63) and linear 
controller (III. 65). The state matrix F of this feedback system has the following form: 

Ax BV c 3 C j BC c i BV c2 

A 2 ^.3 0 0 

B c i + B c 3 C\ B c 3 C 2 A c i B c 2 

C\ C 2 0 0 

Next, we linearize the feedback interconnection of the plant Q and the con- 
troller C . However, in order to do that, first we must determine the values of the 



F : = 



(III. 66) 
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controller states X c \ and X c2 along the trajectory Pc 6 S. From equation III. 62, we 
obtain: 



jX cU = AMX'U + BMljVj j»l (Sic - ec'A C ) T ] T 

+BME + B c3 (q,)j t E, 

~j~X C 2 c = C c \(q c )X c \ c + T> C 2 (q c )E, 

dt 

U c = X c2c . (III.67) 

Notice, since along Pc € £: 



we get 



E 


= 0, 


V c 


= constant, 


flc 


— constant, 


Qc Ac 


= 0, 


X c2c 


= f/ c , a constant, 


d x 

dt cU 


*^-cl^Clc5 


0 


’ — CdXclc- 



(III. 68) 



Now, using Assumption A3, we conclude that 






Cl c 



= 0. 



(III. 69) 



In order to compute the linearization of the feedback interconnection of G 
and C (T(Q ,C )) along Pc 6 £, first observe that J-(Q ,C ) is equal to the feedback 
interconnection of Ge and the system consisting of X(q) and the integrator X&. The 
linearization of Ge along Pc G £ is given by (III. 59). The linearization of the system 
consisting of X(q) and X c2 can be obtained using the trimming values of and X C 2- 
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Define 



U == 

dt 





zu : = 



iT 



dt^E M 



*5 Af 2 £ T |£ t 9c 



nT 



then, along Pc 

-o = 0 1/J 0 0 0 <? Co 

and the linearization of (III. 67) can be written as 

d tv _ 3(Ai(fc)*«i.) , , 8(Boi(fc)£#, , d(B*(q c )E) 

7t cU ~ 55 U+ ~ U + 



(111.70) 



dr? 



3g7 



^0 



d(S a (qc)£E) 



dm 



\VJq ) 



d, v _ 3(C c , too)*,.,), , 3(IMfc)£) 

5i 4 - Yc2 ‘ _ 55 u + 






|ra 0 



Wc = -g^-lw 



(III. 71) 



Consider the expansion of the first term in (III. 71), 
d(A cl (q c )X clc ). 



dm 



l^o 



d(A cl (q c )X clc ) l d(A cl (q c )X cl ) 



dX cl 
d(Ac i(g e )X elc ) 
+ dE 
— A cl ( qco )8X cU . 



\^0 



l^o + 



dX c2 

a(Ai(<? e )x cl ,) 

dE 



Ro + 



d(A cl (g c )X cle ) 

d(f t 0) 



Ro 



d(A cl (q c )X cU ) 

Ro H o |w 0 ) 



dq c 



(III. 72) 



Since X c \ c is zero along the trajectory, their are no extra terms due to the scheduling 
parameter, q c . Similar results are obtained for the terms B c i(q c ), C c i(q c ), and V c2 (q c )- 
Noting that 



R T = 

dt 



Jt^Vg ^8Q E ^8 A e 



iT 



the linearization of (111.67) has the following form: 



j8X clc = A cl 6X clc + B ci dsvi U<d T E ±8A T E ) T + B c2 8E + B c3 ^8E 



dt 

d 



l dt 



dt 



dt 



'dt 



—SX c2 — C c i6X c i c + V c2 8E, 
8U = 8X c2c . 



(III. 73) 
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The state matrix M of Ti(Q , C ) is calculated as 




The proof of the first part of the theorem can now be shown by following the proof 
of Theorem 4.1 in [Ref. 27]. From Assumption A2, it follows that the matrix 

■Acl &c2 

B C 1 T>c2 

is invertible. Therefore, we will use this to define a very useful similarity transforma- 
tion to show that F and M have the same the eigenvalues. Let 



(III. 74) 



-d c l 


B c2 


-1 


X 


Y 


Bci 


V c2 _ 




z 


w 



and set 



T := 



-X 



-Z 



I 




0 0 


r 
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Bel B c 3 
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(III. 75) 



Expression III. 74 implies 



A. C \X + B c 2 
Cd Y + T> C 2 
AciY + B c 2 
Cd A' + T>c2 



/, 

/, 

0 , 

0. 



(111.76) 
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A simple computation using (III. 76) and (III. 75) verifies that 



7 1— 1 



Bd B* 



I 0 
C\ Ci 



0 



0 0 
X Y 
Z IV 



(III. 77) 



and that F — T AIT -1 . Thus, F and M have the same the eigenvalues. 

In order to prove the second part of the theorem, it will suffice to show that 
the linear controllers (III. 73) and (III. 60) have the same transfer function from 



E = 



8v 8y 8z 



8v c 8y c 8z c 



and 



0 = 



8V t 8Q t 8A t 



t T 



(III. 78) 



to 8U . A direct calculation of the Laplace transforms show that 



U(s) = C c i(sl - A c i)~' 
B c2 E(s) 



8V T (s ) 8Q T (s) 8A T (s) 



+ 



5 

T> c 2 E(s) 



+ BczE{s) j + T> c \ 
+ T>czE(s), 



8V T (s) 8n T (s ) 8A t (s ) 



(III. 79) 



for both controllers. ■ 

Thus, the eigenvalues of the linearizations along each trajectory in £ are pre- 
served. Furthermore, the input-output behavior of the linearized operators is pre- 
served in a well-defined sense. The reader is referred to [Ref. 27] for a complete 
discussion on approximations to this method that avoid using pure differentiation. 



3. Implementation Procedure 

The Theorem D.l can be used as follows: first, determine the dynamics of 
the vehicle Q and the set of trimming trajectories £ the vehicle is required to track. 
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This set is parameterized by the range of trimming velocities (u c ), desired flight 
path angles ( 7 c), and desired heading rates (0 C ). Recall, these variable constitute 
a gain-scheduling vector q c = [v c 7 C < 0 C ] T . Next, rewrite the vehicle dynamics 
using generalized error coordinates V#, H#, A #, and Pe to obtain the vehicle’s error 
dynamics Qe- Linearize Qe about a finite number of trajectories in E. Use these linear 
models to design a finite number ( k ) of trajectory tracking controllers {Ci i = 1, k}. 
Gain schedule = 1 , A:} utilizing a favorite interpolation or gain-scheduling 

technique to obtain the linear gain-scheduled controller C\ (g c ). Now, implement this 
controller on the nonlinear plant Q according to the expression III.62. 

It is worth emphasizing the following important properties of the controller C : 

• The result in Theorem D.l holds for all trajectories in E. 

• The structure of the controller C is easily obtained from that of the linear 
controllers. 

• Since all the closed-loop transfer functions of the local linearizations are pre- 
served, at the level of local linear analysis, the controller does not introduce 
any additional noise amplification despite the presence of a differentiation op- 
erator. 

• Along trajectories Pc G £, X C 2 c = U c and X clc = 0. Therefore, the trimming 
values of the control inputs are naturally provided by the integrator block with 
state X C 2 c - 

• The integrators X C 2 are directly at the input of the plant, which makes it 
straightforward to implement anti-windup schemes. This becomes necessary 
in applications where the input U is hard limited due to actuator saturation, 
for example. 

• The inputs to the controller, Pc and A^, can be computed directly from the 
vector q c . 

• The trim values and U c are not required in the controller implemen- 

tation. 

• Along the trajectories Pc G £, the controller guarantees that the steady state 
value of error vector E is zero, which follows from the fact that the controller 
solves the controller implementation problem. This is in sharp contrast to 
standard LOS guidance schemes. 
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E. EXAMPLE 



In this section we apply the methodology developed in Section D to the design 
and implementation of an integrated guidance and control system for a fixed wing 
unmanned air vehicle. The vehicle’s stability and control derivatives are the subject 
of a section of chapter IV. Using the notation in Section B, the Frog dynamics have 
the following form: 



G = 



±v =^v(v,n,A) + Mv,n)H(v,n,u), 



= Fa(V,n,A) + la(V,n)H(V t Sl,U) t 
i p = ! bKV, 



(III. 80) 



i; A =QSl. 



Following the development in Section B, the set of trimming trajectories £ for 
the vehicle is defined as follows: 



i p c=b nv c , 

21 Ac = Qc ^c, 

Fv(V c ,n c ,Ac) + lv(VcMH(Vc,Slc,U e ) = 0, 

Fn{Vc, &c, Ac) + Za(Vc, flc)W(Vc, £lc, U c ) — 0 , J 

(III. 81) 

where Pc and Ac can be computed using the scheduling vector q — [u c 7 c ip c ] T . 
Now, given [v c 7 c i/> c ] T and (3 C = 0, we can solve for Vc,Hc, U c , and Ac: 



’ Pc 
Ac 



Tv(v c ,n c ,A c) + iv(v c ,n c )H(Vc,nc,u c ) = o, 

Ta(V c ,n c ,Ac) + MVc,nc)H(V Ci nc,U c ) = 0, 

Qc'^c — A c — 0 
IIVoll - », = 0, 

A - sin- 1 £ = 0, 

Vt 

7 c — [0 1 0]arg(^nfn) = 0, (III.82) 

where the arg function extracts the angles X from the rotation matrix 1Z{X): X = 
arg(K{X)). 
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Using the solution to equations III. 82, the linear model for the vehicle rep- 
resented by equations III. 59 was obtained along a trajectory characterized by the 
velocity of 73 feet per second, flight path angle of zero, and heading rate of 10 degrees 
per second. This model was used to design a linear trajectory tracking controller for 
the vehicle. 



1. Design Requirements and Linear Controller De- 
sign 

Design requirements for the linear trajectory tracking controller included 

• Zero Steady State Error: Achieve zero steady state tracking errors of all 
trajectories in £. Achieve zero steady state tracking error of indicated airspeed 
while on any trajectory in £. 

• Bandwidth Requirements: The command-loop bandwidth for each com- 
mand channel should be no greater than 1 radian per second and no less than 
1/10 radian per second; the control-loop bandwidth should not exceed 12 ra- 
dians per second for the elevator, aileron and rudder loops, and 5 radians per 
second for the throttle loop. These numbers represent 50% of the correspond- 
ing actuator bandwidths, and shall ensure that the actuators are not driven 
beyond their linear operating range. 

• Closed Loop Damping and Stability Margins: The dominant closed- 
loop eigenvalues should have a damping ratio of at least 0.5. Simultaneous 
gain and phase margins of 6db and 45 degrees in each control loop must be 
achieved. 

The methodology selected for linear control system design was synthe- 

sis [Ref. 12]. This method rests on a firm theoretical basis, and leads naturally to 
an interpretation of control design specifications in the frequency domain. Further- 
more, it provides clear guidelines for the design of controllers so as to achieve robust 
performance in the presence of plant uncertainty. The basic steps in the controller- 
design procedure, including the development of the synthesis model, were done us- 
ing the approach described in [Ref. 24]. This approach provides an intuitive and 
straightforward way for converting the design requirements into the weights for the 
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'H 00 synthesis model. Consider Figure 35. Here C; is the controller to be designed, 
Qi is the linear model of the vehicle. 




Figure 35. Synthesis Model 



In Figure 35, the vector of exogenous inputs w represents the commanded 
inputs. The vector y\ represents lateral and vertical displacement states of the linear 
model as well as the vehicle’s velocity. The regulated output 2 includes the outputs 
of the weighting matrices W\ and W 2 . These matrices had the following form: 



W x 



W 2 



^00 

5 

0 f 0 , 

0 0 a 

5 

c 4 0 0 

0 C5 0 , 

0 0 c 6 



(III. 83) 



where the constants c 4 -, i = 1,6 were used as the design knobs adjusted to meet the 
closed-loop tracking, damping, control, and command loop bandwidth requirements. 
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Values of W\ and W 2 were iterated on, and used to obtain a linear trajectory 
tracking controller: 



8E = [ 8v 8y 8z] T — [ 8v c 8y c 8z c ] T , 



^8X C 1 — A c \ 8X c i + B C \8X C 2 , 

8U = q f{C cl 8X cl + V cl [8V T 8tt T 8 A t e ] t + V c2 8X c2 + V* 8E), 
where the TYoo state-feedback gain is X = [ C c \ X > C 2 The feedback system 

consisting of the plant Q\ and the controller Ci was found to meet all the design 
specifications given earlier in this section. Since the control surface effectiveness is 
proportional to the dynamic pressure, the controller C\ was gain-scheduled on q. 
The value <jo represents the nominal value of q. 



Cl (q , ) 



2. Implementation and Simulation Results 



Using the formulae provided in Section D, the family of linear gain-scheduled 
controllers Ci (q) was implemented on the nonlinear plant Q as follows: 



C := 



[0 y z\ 


=f n(p - Pc(s 0 )), 


E 


= [v-v c y z]. 


X c \ 


= Ac lX c i + B C 2 E , 


X c2 


= f{C cl X cl +V cl {±V T 


u 


= x c2 . 



dt w 



} dt l 






We emphasize that the implementation equations for the controller C do not require 
the computation of flc, U c and Vc- Moreover, since c = [0 0 ipc] 7 ', the controller 
must only be provided with ip c and Pq when steering the aircraft along the trajectory. 
These are the critical advantages of the proposed methodology. The acceleration 
term can be computed using onboard sensors without resorting to differentiation. 
Therefore, the only term which could not be computed directly was In this case, 
the differentiation operator ^ was replaced by a causal operator with the transfer 
function [Ref. 27]. 
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The trajectory tracking controller developed above was tested using a number 
of trajectories. One such trajectory consisted of a straight line transitioning into a 
helix shown in Figure 36. This trajectory is characterized by a cruise velocity of 73 
feet per second. Initially, the trajectory is aligned with the inertial x-axis. After 
proceeding along the x axis for 3000 feet, the trajectory turns into a helix with a 
radius of 1000 feet and climb angle of 5 degrees. Consider Figure 37, which shows 
the time history of the position error, bank and pitch angles, and indicated airspeed 
along the trajectory. Clearly, the controller drives the vehicle along this trajectory 
with zero steady state position errors while maintaining 73 fps indicated airspeed. 




Figure 36. The trajectory tracked in simulation. 



F. CONCLUSIONS 

A new method was introduced for designing and implementing integrated guid- 
ance and control systems for autonomous vehicles. The starting point is a family of 
linear controllers with integral action designed for linearizations of the nonlinear equa- 
tions of motion described in an appropriate state space. Based on this family, the 
method produces a gain scheduled controller that preserves the input-output proper- 
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Figure 37. Time history of position errors, Euler angles and airspeed along the 
trajectory. 



ties of the original linear closed-loop systems as well as the closed-loop eigenvalues. 
The key feature of the method is the ability to automatically reconfigure the control 
inputs of the vehicle to provide for proper control action as the body tracks an inertial 
trajectory in free space while maintaining constant airspeed. The method is simple 
to apply and leads to a nonlinear controller with a structure similar to that of the 
original linear design. 
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IV. 



RAPID PROTOTYPING SYSTEM FOR 
FLIGHT TEST OF AN UAV 

A. INTRODUCTION 

This chapter describes the development of a rapid prototyping system for flight 
testing of guidance, navigation, and control algorithms for unmanned air vehicles. The 
system affords a small team the ability to take a new concept in guidance, navigation, 
and control from initial conception to flight test. In order to do this, a number 
of engineering problems had to be overcome addressing a gamut of issues including 
weight, power, portability, risk, electronic interference, vibration, manpower, etc. The 
main contribution of the project is the proof of concept flight test demonstration of a 
new integrated guidance and control algorithm (see Chapter III). The success of this 
endeavor had a synergistic effect on several other projects, paving the way for joint 
ventures in voice controlled flight, coastal mapping, and autonomous landing. The 
project is viewed as the foundation of a long-term investment in innovative unmanned 
air vehicle applications leveraging, in part, the operational experience of the officer 
students in the avionics curriculum. 

Testing of a new algorithm, sensor package, vehicle, etc., requires expertise 
from many branches of the engineering sciences, especially aeronautic, electrical, and 
computer. It is potentially costly and time consuming, as well as having the potential 
for catastrophic failure. When successfully done, however, it provides developmental 
information, insight and data that are unavailable from other sources. All of our 
theoretical and numerical results must be verified by some form of experiment, and 
flight testing is often the best way to do this. 

The chapter begins with a conceptual discussion of the Rapid Flight Test Pro- 
totyping System (RFTPS). Motivation for its development is addressed. Next, a 
description of the hardware components is given. The main contribution of the chap- 
ter is discussed in the last section, where an application of the RFTPS to the problem 
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of integrated guidance and control introduced in Chapter III is presented. The full 
capabilities of the RFTPS are demonstrated when this novel guidance algorithm is 
taken from theoretical development to flight test. 

B. SYSTEM DESCRIPTION 

The RFTPS consists of a test bed unmanned air vehicle equipped with a com- 
plete avionics suite necessary for an autonomous flight as shown in Figure 38, and of 
a ground station responsible for flight control of the UAV and flight data collection as 
shown in Figure 39. A functional block diagram of the RFTPS is shown in Figure 40. 
The key goal was to use off-the shelf-technology as much as possible, thus exploiting 
the economy of scale of a number of commercial industries. Furthermore, if the UAV 
development program is to span many years, and to draw on the talents of the officer 
students in the future, the RFTPS had to emphasize high level algorithm design. Low 
level code and device driver generation is kept to a minimum with the vast majority 
of the code "writing” being done via autocode tools. The system architecture is open, 
providing the ability to add, remove or change real time input/output (I/O). Compu- 
tational power can be increased as mission requirements dictate. The telemetry links 
are secure, yet low power and unobtrusive to the public, not requiring advance per- 
mission for use, special frequencies from a government authority, or special airspace. 
The onboard components are light weight and low power, allowing for the inclusion 
of additional payload. 

1. RFTPS Capabilities 

The RFTPS developed provides the following capabilities. 

• Within the RFTPS environment, one can synthesize, analyze and simulate 
guidance, navigation, control, and mission management algorithms using a 
high level development language. The same code that ran the simulation, flies 
the vehicle. 

• Algorithms are seamlessly moved from the high level design and simulation 
environment to the real time processor. 
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Figure 38. Frog : The unmanned air vehicle at the Naval Postgraduate School. 



• The RFTPS utilizes industry standard I/O including digital to analog, analog 
to digital, serial, and pulse width modulation capabilities. 

• The RFTPS is portable, easily fitting in a car. In general, testing will occur 
at fields away from the immediate vicinity of the Naval Postgraduate School. 

• The unmanned air vehicle can be flown manually, autonomously, or using a 
combination of the two. For instance, automatic control of the lateral axis can 
be tested while the elevator and throttle are controlled manually. 

• All I/O and internal algorithm variables can be monitored, collected and an- 
alyzed within the RFTPS environment. 

2. Cost, Safety and Other Considerations 

Cost and risk are two leading, and at times competing, concerns that had to 
be effectively handled. Since initial testing is to occur within line of sight at all times, 
a pulse width modulated (PWM) remote control system manufactured by Futaba was 
chosen. Testing of a new control algorithm is similar to handing over control of the 
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Figure 39. The base station of the RFTPS in use at the airfield. 



aircraft to a student pilot. The algorithm should have full freedom to perform, yet 
adequate safeguards must exist in case it fails. With some modifications, the extensive 
master-slave flight training capabilities built in to the existing RC transmitters were 
exploited. A significant portion of the cost of the RFTPS resides in the real time 
processor, I/O board and modules, and in the host computer. Additionally, while 
compact, the weight and power requirements of these components are significant 
when compared to onboard power and payload available. In order to gain additional 
payload, and in order to manage the risk associated with the loss of an expensive 
computer package, the real time controller was kept on the ground. Sensor and 
control links to the real time controller were bridged via RF components described 
later. 
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Figure 40. RFTPS Hardware Architecture 



3. Components 

The centerpiece of the RFTPS ground station is the AC100/C30 system from 
Integrated Systems Incorporated. The key feature of this product is its autocode 
tools. With a relatively short time available for research by the officer students, em- 
phasis had to be shifted from code writing, debugging and maintenance to algorithm 
development. AC100/C30 utilizes ’’Xmath/SystemBuild”, a graphical programming 
environment that uses a high level block diagram paradigm for modeling of linear 
and nonlinear systems. Within the ’’Xmath/SystemBuild” environment, the algo- 
rithm can be built, simulated, tested, and debugged. Real-time code can then be 
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generated for execution on the real time processor. 

Currently, the ”Xmath/SystemBuild” environment resides on a Sun worksta- 
tion. PC’s running Windows NT are also supported. This allows for a significant 
downsizing of the ground station should resources become available in the future. 
Communication with the real time processor is via an Ethernet bus using TCP/IP 
protocol. AC100/C30 provides excellent animation tools for building graphical user 
interfaces (GUI). Through the appropriate design of these interfaces, the flight test 
team can monitor, modify ,and control the actions of the real time processor. The 
GUI resides on the workstation. Communication between the workstation and the 
real time processor is via the Ethernet connection, and is managed by a host PC. Ad- 
ditionally, the host PC provides power to the real time processor, as well as providing 
utilities for compiling, linking, and downloading the C-code. 

The I/O consists of four multi-mode, bi-directional, serial ports utilizing RS- 
232 protocol, a 16 channel pulse width modulation port capable of measuring up to 
sixteen PWM signals or generating up to six PWM signals, and a six channel digital- 
to- analog converter. The I/O modules are hosted by the same PC that holds the 
real time processor. The real time processor is a single Texas Instruments Digital 
Signal Processor (TMS320C30). The capability exists for upgrading the processor, 
or running multiple processors, to meet computational demands of future projects. 

The control configuration of the air vehicle is conventional with three indepen- 
dent surfaces (elevator, aileron, rudder) and a throttle. Manual control is provided 
via a Futaba, dual conversion, PWM transmitter utilizing the portion of the radio 
spectrum reserved for Radio Controlled (RC) flight, 72.030 MHz to 72.990 MHz. 
Precautions entail a search of the electronic spectrum utilizing a hand held spectrum 
analyzer, as well as standard procedures employed by RC hobbyists to avoid two in- 
dividuals selecting the same frequency locally. Built in capabilities of the transmitter 
include the ability to transmit one or more signals from a slave transmitter. The 
slave transmitter is a modified Futaba transmitter where the manual control effectors 
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have been replaced by a direct connection to the digital-to-analog I/O module. In 
this way, an exogenous source (RFTPS) can be given control of one, some, or all of 
the control actuators of the aircraft using the same RC link currently controlling the 
aircraft. 

The sensor suite onboard the air vehicle consists of an inertial measurement 
unit (IMU) with a three axis rate gyro, three axis accelerometer, magnetic heading 
indicator, two axis pendulum, phase differential GPS receiver, elevator, aileron, and 
rudder actuator position sensors, and angle of attack, side-slip angle, pitot-static, and 
static pressure air data sensors. A four channel analog-to-digital converter is used to 
capture any four of the following sensors: elevator, aileron, rudder actuator position, 
angle of attack, side-slip angle, dynamic pressure, or static pressure. 

Communication between the sensors and the real time processor is via a serial 
link. Low power, matched, spread spectrum RF links provide up to 115 Kbaud rates 
at over 10 miles range. They require no license, and can be used anywhere in the 
United States. Additionally, the onboard GPS unit maintains contact via a serial link 
with a GPS receiver on the ground, which provides differential corrections. 

Power for the onboard avionics is supplied from a lithium-ion battery. The 
battery provides 6 hours of continuous use. The power budget of the onboard com- 
ponents is shown in Table III. 





Voltage 


Current 


Power 


IMU 


24 Volts 


400 Milliamps 


9.6 Milliwatts 


GPS 


12 Volts 


250 Milliamps 


3.0 Milliwatts 


Telemetry 1 


24 Volts 


300 Milliamps 


7.2 Milliwatts 


Telemetry 2 


24 Volts 


300 Milliamps 


7.2 Milliwatts 


Total Req. 






27 Milliwatts 


Battery 






> 200 mW-Hrs 



Table III. Power Budget of Onboard Components 
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C. TRAJECTORY CONTROL: AN APPLICATION 

The intent of this application was to demonstrate the utility of the RFTPS by 
flight testing a new integrated guidance and control algorithm that was shown to have 
good performance and robustness properties both theoretically and in simulation (see 
II). This project was chosen because, on one hand, it is a totally unique application 
in terms of the guidance and control algorithms implemented. On the other hand, 
to implement this algorithm, most of the steps required of any application involving 
autonomous flight of an air vehicle must be accomplished. If done correctly, the 
foundation will be laid for follow on projects and joint ventures utilizing unmanned 
air vehicles. 

Generic tasks fundamental to guidance, navigation and control algorithm de- 
velopment were accomplished first. To begin with, a high fidelity model of the test 
bed vehicle, nicknamed Frog , was developed. This involved a complete lateral and 
longitudinal parameter identification of Frog’s stability and control derivatives. This 
lead to the development of a six degree of freedom, nonlinear simulation. In order 
to manage the risk associated with the autonomous flight, it was decided to add an 
onboard inner-loop controller. The inner-loop controller modifies the vehicle dynam- 
ics such that displacements of the primary longitudinal and lateral control effectors 
correspond to various trimming trajectories , as defined in Chapter III., The inner-loop 
controller selected was a commercial autopilot whose control laws had to be identified 
as well. Additionally, accurate and timely calibration of the I/O needed to become 
an integral part of the RFTPS to account for changing environmental conditions, 
battery levels, etc. 

This section begins with a brief background of the parameter estimation algo- 
rithm used to identify a model of Frog. Next, the experimental testing completed to 
obtain data used in the parameter identification process is summarized. Flight test 
data is compared with the simulation results to demonstrate the efficacy of the ef- 
fort. Then, a linear controller design based on modeling results obtained is discussed. 
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Next, a nonlinear controller implemention based on the theory developed in Chapter 
III is presented. Penally, flight tests conducted using the RFTPS to track an inter- 
esting representative trajectory are discussed. Data from flight tests are presented 
and contrasted with results from simulations. Along the way, implementation issues 
germane to the flight testing process are explained. 

1. Theoretical Background for Model Identification 

Parameter identification of the test bed vehicle and inner-loop autopilot was 
accomplished based on a series of test flights and lab experiments. The Maximum 
Likelihood Method was chosen to identify the parameters of the model. The back- 
ground is extracted primarily from [Ref. 18]. Assume that the model to be identified 
has the following form 

x = A(p)x + B(p)u, 
y = C(p)x + D(p)u, 
x = state vector, 
u = input vector, 
y — otput vector, 
p = parameter vector. 

The problem at hand is to estimate the parameter vector, p, given a data set of 
measured inputs, u, and outputs, y. To that end, denote the estimated output based 
on the parameter vector p as p(p), and use it to define a cost function 

J(p) = II y(p)-y\\l 

The method seeks to minimize the cost, J, by varying the parameter vector p. Given 

an initial guess of the parameter vector p 0 , the Jacobean of the cost function J is 

calculated by numerically perturbing the parameter vector p 0 

dJ_ _ yjpo + Sp) - yjpo ) 
dp 6p ’ 

=: Ho. 
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Then, a linear approximation for y(p + 8p) is 

y(p + Gp) ~ v(Po) + H 0 8p. 

The cost J can now be expressed as a function of the perturbation of the parameter 
vector as 

J(Sp) = \\«(p+6p)-«\\l (IV.l) 

= l|y(Po) + HoSp — y\\l, 

= IIjKpo) - y 111 + 2y(p 0 ) T H 0 6p - 28p T Hj : y + Sp T H„8p. 

In order to determine a 8p which minimizes equation IV.l, set 

dJ_ 

d(8p) 

from which it is evident that 

HlH 0 6p = HZ(y-y(p 0 )) (IV.2) 

must hold. 

A common problem encountered in solving a parameter identification problem 
is the presence of redundant parameters in the model. This results in the matrix 
Hq H 0 being singular, and presents a problem if a solution to equation IV.2 is sought 
in terms of 8p. In order to avoid this singularity, the following problem is solved 
instead, 

(HlH 0 + pI)8p = HZ(y-y(p 0 )). (IV.3) 

With the descent direction 8p known, a line search is used to calculate the step size. 
This method is guaranteed to converge to a minimum for the cost function J, although 
not necessarily a global minimum [Ref. 18]. 

The appealing feature of this method for the problem at hand is that it can 
exploit a great deal information known about the model. For instance, the order of the 



2 H% y( Po ) - 2 Hi y + 2// 0 T H 0 8p, 
0 , 
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plant and most of the structure in the system matrix is known, as is a reasonable range 
of values for the parameter vector p based on empirical methods (USAF DATCOM). 

2. Vehicle Model Identification 

The conventional configuration of Frog suggested that the parameter identifi- 
cation problem could be decoupled into longitudinal and lateral dynamics. Maximum 
likelihood parameter identification was used to refine existing analytic estimates of 
the longitudinal stability and control derivatives [Ref. 34]. The results for a few key 
longitudinal stability derivatives are compared to analytic estimates in Table IV. The 
primary difference was in the estimate of the elevator effectiveness, which was less 
than half as effective has previously thought. 



Derivative 


Analytic 


Experimental 


%A 


Cl. 


4.30 


4.09 


5.13 


Cm q 


-.417 


-.557 


-25.1 


C L 6e 


-1.12 


-.391 


186 


C M (e 


-1.62 


-1.05 


54.3 



Table IV. Comparison of selected longitudinal derivatives. 



Figure 41 shows the response, in simulation, of a longitudinal model of Frog 
using these derivatives to an elevator doublet. Comparisons are made to the results 
from flight test. 

Next, the process was repeated for the lateral axis. Aileron and rudder dou- 
blets were executed while aileron and rudder position, roll and yaw rates, and side-slip 
angle were measured. Results, for a few key lateral stability derivatives, are compared 
to analytic estimates in Table V. The biggest difference was the value of which 
increased by over 300 percent. This doubled the estimate of the natural frequency of 
the dutch roll mode. Figure 42 shows the response, in simulation, of a lateral model 
of Frog using these derivatives to an aileron doublet. Comparison to the results from 
flight test are also shown. 
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Figure 41. An elevator doublet is used to excite the longitudinal dynamics of Frog. 
Test flight data is compared with simulation results to assess validity of the model. 



3. Autopilot Model Identification 

The inner-loop autopilot is a “black box” containing both sensors and con- 
troller logic. The intent was to model the unit as closely as possible without disas- 
sembling it. The published function of the autopilot is to control vertical speed and 
turn rate of the aircraft. This naturally decouples into an identification problem for 
the lateral channel, and an identification problem for the longitudinal channel. 

The lateral channel employs a rate gyro to track turn rate commands via 
feedback to the aileron. The general structure is shown in Figure 43. In order to 
identify the dynamics of the block labeled T\ in Figure 43, the unit was rotated at 
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Derivative 


Analytic 


Experimental 


%A 


cy* 


-.310 


-.987 


-68.6 


Ci B 


-.051 


-.094 


-45.7 


Cn 3 


.058 


.176 


-67.0 


Ci la 


.181 


.239 


-24.3 



Table V. Comparison of selected lateral derivatives. 



differing yaw rates. This provided a variable input signal to the feedback path with a 
frequency content covering 0 to 20 radians per second. The commanded yaw rate, r c , 
was held constant. The feedback loop was broken at the summing junction of r c and 
r/6, and the feedback signal was captured. Parameter identification algorithms were 
used to determine that the feedback loop could be approximated by the following 
transfer function: 



T, 



~.ls + 1 
5 + 1 



(IV.4) 



With the autopilot on, a step command in turn rate was transmitted to the 
vehicle in flight. Vehicle turn rate, as measured by the IMU, was recorded and used 
to estimate the value of Ki at at 0.25. Test data from that flight is shown in Figure 44 
and compared with simulation results using the autopilot model. 

The longitudinal channel of the autopilot senses the rate of change of static 
pressure in order to control the vehicle’s vertical velocity via feedback to the elevator. 
Flight tests data capturing the response of the vehicle to a step input in climb rate 
command was used for model identification. The identified model of the longitudinal 
channel of the autopilot is shown in Figure 45. 

With accurate models of Frog and of the onboard autopilot identified, it was 
a simple matter to determine the command bandwidths of the inner-loop controller 
(see Figure 46). For this project, this would be the limiting factor in achievable 
performance of the guidance algorithm. Hardware-in-the-loop testing of the actua- 
tors found their bandwidth to be an order of magnitude higher than the inner-loop 
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time (seconds) 



Figure 42. An aileron doublet is used to excite the lateral dynamics of Frog. Test 
flight data is compared with simulation results to asses validity of the model. 

command bandwidths. A natural next step would be to remove the autopilot after 
sufficient time and experience have reduced the risk exposure to acceptable levels, 
and exploit the extra actuator bandwidth available. 

4. Design of the Controller 

The design requirements to be met were as follows: 

1. Tracking Requirements 

• The controller should achieve perfect steady state tracking of trajectories 
in E in the presence of a constant disturbance. For a definition of the set 
£, see Chapter III. 
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Figure 43. Block diagram of the lateral channel of the inner-loop autopilot. 



2. Bandwidth Requirements 

• Command bandwidths along the lateral and longitudinal channel should be 
maximized to ensure tight tracking of the trajectory. Adequate frequency 
separation between the inner and outer loops should be assured. 

3. Closed Loop Damping and Stability Margins 




limo (saconde) 



Figure 44. Flight test data is used to identify the dynamics of the autopilot. A step 
signal is sent to the lateral channel of the autopilot. Measured yaw rate is compared 
to simulation results. 
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Figure 45. Block diagram of the longitudinal channel of the inner-loop autopilot. 



• The dominant closed-loop eigenvalues should have a damping ration of at 
least 0.5. Phase and gain margins should be no less than 45 degrees and 6 
dB respectively in each control loop. 

4. Implementation Requirements 

• Control of Frog from the console, open-loop, must be possible in order to 
maneuver the vehicle to a suitable location overhead the field based on 
visual cues and displayed navigational data. 




Frequency [rad/sec] 



Figure 46. Bode plot of the yaw rate command to yaw rate and climb rate command 
to climb rate for Frog with the autopilot on. 
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• Switching from open-loop control to autonomous flight should be accom- 
plished from the console at the discretion of the console operator. The 
switch should be bumpless with no undesirable transients. 

• The definition of the inertial trajectory passed to the controller should be 
easily defined from the console, and expressed in terms of the parameters 
introduced in Chapter III, namely, helix angle ( 7 C ), radius (6 C ), and turn 
rate 



a. Linear Controller Synthesis 

The synthesis of the controller was an application of the theory outlined 
in Chapter III where the trimming trajectories were parameterized by the vector 
7 Jc = [v c ipc 7 c] 6 7£ 3 . This was useful for the derivation of the generalized error 

dynamics. For the application at hand, however, it was convenient to remove the 
explicit dependence on v c . Since along Pc € £ 

_ t; c cos( 7 c ) 



there is no problem in using fj c = b c t/> c 7 c j € TZ 3 to parameterize Pc G £ vice r/ c . 
This was done strictly for convenience, since during the flight test, the throttle was 
to be left at a cruise power setting, and the airspeed was not explicitly controlled. 
Flight tests tracked a single trajectory in £ identified as 

1 T 

fjc = 1146 5 0 , (IV. 5) 



where the units are feet, degrees and seconds, respectively. Based on (III. 6), (III. 7) 
and the model identification results in section 2 and 3, the nonlinear plant, 

' Tt v =*V(v,n,A) + Tv(v,n)«(v;n,t/), 

g = I =/h(v,n,A) + Mv,n)H(v,n,u), 

Tt P =bRV, 

, it A 

was available for controller synthesis. Given (IV. 5) and (IV. 6), the trimming values, 
Vc, Dc, and Uc were solved as the solution to 



fv(Vc,n c Ac) + iv(Vc,nc)H(Vc,n c ,u c ) = o, 
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Xn(Vc, Clc, Ac) + Zq{Vc, VtcYHiVc, &c, U c ) = 0, 

<2 _1 i7c — Ac = 0, 



11^11 — 



^cV’c 

cos( 7 c) 



= 0 



fie - sin 1 



bc^c 

cos(7c) 

Vt 



= 0 , 



7c -[0 1 0 ] arg^n'eR) = 0 , 



(IV.7) 



where the arg function extracts the angles X from the rotation matrix 1Z(X): X = 
arg(lZ(X)). Numerically solving (IV.7), the trim values were found to be 



Vn = 



n T c = 



A l c = 



87.7 2.037 2.43 
0.00176 0.029 0.087 

- \T 

13.41 0.023 0 



where, as before, standard units are feet, degrees and seconds. 

Based on lab and flight testing of the onboard sensors, the decision was 
made to design an output feedback controller of the form (see Figure IV. 8): 



8E = [Sy Sz] T - [8y c 8z c ] T , 



%8X cl = A cl 8X cl +B c2 8X c2 + B c3 8E, 

&X* = 6E, 

8U = C c \ £A c i + T) C 28 X C 2 + 'D&8E, 



(IV.8) 



where 8U = [r c h c . 

Classical control techniques were employed to design the controller, C \ , 
shown in Figure 47. Performance objectives included setting the loop gain crossover 
frequencies at 0.5 radians per second along each channel based on the inner-loop 
command bandwidths (see Figure 46). Root-locus analysis along the lateral channel 
provides the most intuitive view of the control problem (see Figure 48). The open- 
loop plant has unstable zeros along the lateral channel due to the dihedral effect of 
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Figure 47. Linear Controller Design 

the high wing design. These unstable zeros will naturally attract the free integrators 
in the synthesis model if left unattended. 




Figure 48. Root-locus for the lateral channel. 

The straightforward solution was to interlace the pole-zero structure 
of the synthesis model with a stable pair of complex zeros sufficiently close to the 
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free integrators. The robustness limit was set by the dutch roll poles which were 
forced to eventually migrate to the non-minimum phase zeros. Modification of the 
dutch roll behavior would require a different approach than that chosen with the 
inclusion of the inner-loop autopilot. Replacement of the inner-loop autopilot is a 
natural extension for future flight test work. Design concerns and solutions along 
the longitudinal channel were qualitatively similar. The resulting design was a sixth 
order controller. 

Nyquist diagrams of the loop transfer function for the lateral and lon- 
gitudinal channels are shown in Figure 49 and 50. There it can be seen that phase 
and gain margin design requirements along each channel have been met. 




Real 



Figure 49. Nyquist diagram of the loop transfer function for the lateral channel. The 
phase margin is 112 degrees and the gain margin is 6 dB. 

5. Nonlinear Implementation 

The controller designed in the preceding section was implemented on the non- 
linear plant using the implementation proposed in Chapter III and is given next (see 
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Figure 50. Nyquist diagram of the loop transfer function for the longitudinal channel. 
The phase margin is 145 degrees and the gain margin is infinite. 



Figure 51): 






[0 y z) T 

E 

d v 
dt Aci 

d y 

Tt A a 

u 



fK(P - Pc(s o)), 

[y-y c z- z c ], 

A cl X cl + B c2 E + Bc 3^E, 
C cl X cl + V c2 E + X>c3 -^E, 

X c2 . 



(IV.9) 



[b c T c V’c] 



P 




Figure 51. Nonlinear Implementation of the Controller 
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The following algorithm was implemented in the block labeled, Calc Me , in 
Figure 51 in order to compute Me- Assume the position of the vehicle is [a:;, yi zi] T , 
and let the center of the trajectory be specified as [x c , y c z c ] T . Then, the projection 
of the position of the vehicle onto the trajectory, as defined by fj c , is 

r i t 



Pc{s o) 



x 0 yo Zq i 



where 



s 0 cos(7 c ) 

x 0 = b c c os( - ) + x c , 

be 

. So cos(7c) 

yo = b c sin( ) + y c , 

be 

z 0 = V> c 6 c s 0 tan(7c) + z c . 



(IV.10) 



When 7 C is zero, the parameter, s 0 , can be computed by considering the projection 
of P onto a co-level circle with radius, 6 C , and center [a: c , y c z{\ T . Then, 

5o = 6 c tan- 1 (^^-). (IV.ll) 

Xi~X c 



The position error resolved in the Frenet frame attached at the point Pc(s o) 



x err 




-cos( 7 c)sin(f^) cos( 7 c )cos(f^) -sin( 7 c ) 




1 

0 

H 

1 


l ferr 


— 


-cos(f^) -sin(f^) 0 




yi - x 0 


Z er r 




sin( 7 c)sin(f^) - sin( 7 c ) cos(f^) cos( 7 c ) 




i 

O 

l 

£ 

\ 



= : M E - (IV. 12) 



If the helix angle is not zero, then a simple solution is not available. The point 
can be found, however, by noting that Me must be of the form [0 y err z err \ . Then 
using (IV. 12), 

- cos( 7 c ) sin(^)(a;/ - x 0 ) + cos( 7 c ) cos(^)(j// - y 0 ) - sm(^ c )(zi - z 0 ) = 0, 
be b c 

(IV. 13) 
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must hold. Using (IV. 10), (IV. 13) is rewritten as 



/ a • , / s o cos(7 c ) 

cos( 7 c ) sin( — )(x/ - b c cos( ) + x c )+ 

be b c 

, x Ow , • r s 0 cos( 7c ). 

cos( 7 c )cos( — )(y, - b c sin( ) + y c )~ 



sin( 7 c )(z/ — V> c 6 c s 0 tan( 7 c ) + z c ) = 0 , 



(IV. 14) 



which can be solved for sq. 



D. CONTROLLER IMPLEMENTATION FOR FLIGHT 
TEST 



During a flight test it is critical that the transition from manual to autonomous 
flight occur without any undesirable transients. The method chosen to achieve this 
was to initiate autonomous flight with the vehicle on the trajectory, and with the 
controller states recruited to their nominal values. An elegant approach toward partial 
completion of this task was to restructure the controller using the ^-Implementation. 
One important convenience of this structure is the straightforward manner in which 
the controller states can be recruited at the initiation of autonomous flight based on 
the trajectory definition parameter fj c . Let t on be the time that the console operator 
turns the controller on. Then, setting 



A c 2(^on) 



A cl (ton) 



V>c cos( 7 c) 
ipcbc tan( 7c ) 

0 

0 

0 

0 



(IV. 15) 



perfectly recruits the controller states at the initiation of autonomous flight. The 
structure in (IV. 9) is also convenient for implementing hard limits on the feedback 
signal. It was experimentally determined that the command signals to the autopilot 
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should be kept within the following ranges 

— 15 degrees per second < b c < +15 degrees per second, 

—2000 feet per minute < h c < +2000 feet per minute. (IV. 16) 

Since the feedback signals are the outputs of integrators in (IV. 9), the states X C 2 were 
easily hard limited to the values in (IV. 16), with anti-windup implemented on the 
associated integrators. 

To affect a smooth transition, all that remained was to initiate autonomous 
flight when the vehicle was on the reference trajectory and aligned with it. Since there 
were no mission requirements regarding the center of the trajectory, the coordinates 
were computed to place the vehicle on the trajectory at initiation of autonomous 
flight. Furthermore, the orientation of the Frenet frame was matched to the vehicle’s 
heading. Let tpon be the vehicle heading at t on . Let the trajectory parameter be 
specified by (IV. 5). Then, setting 





cos(^) cos (7 C ) 


X C X QJl 


be 


1 




sin (xfti) cos 


(7c) 


Vc — Von 


b c 




= Zon 







defines a trajectory where Me at t on , as defined by (IV. 12), is zero, and A 3 = A c 3 , 
where the subscript indicates the third element of the vector. 

1. Additional Implementation Issues 

Prior to autonomous flight of Frog, it was decided to test the RFTPS by flying 
Frog remotely from the workstation console. The console operator was provided with 
the appropriate displays showing Frog s position, heading, and velocity, and used that 
information to command the vehicle’s turn rate and climb/descent rate in order to 
keep Frog in the local operating area. 

The vehicle’s position, as reported by the onboard GPS, is in units of latitude, 
longitude and height above mean sea level. In order to provide a more intuitive 
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navigational reference frame, the location of the workstation at the field where the 
flight tests were conducted was surveyed. This postion was used to define the origin 
of a local tangent plane. Let [Ai, A 2 , h] T be the vehicle’s position, as reported by the 
GPS receiver, where 

Ai = degrees of latitude, 

A 2 = degrees of longitude, 

h = height GPS in meters. (IV. 17) 

The variable h is a height referenced to the surface of an ellipsoid approxi- 
mating the shape of the earth. Two important parameters describing this reference 
ellipsoid, termed WGS8 4, are its eccentricity factor (e) and the length of its semi- 
major axis (a). The distance from the origin of the WGS84 ellipsoid to a point on 
its surface where the local latitude is Ai is given by 

N = T- = ■ =. (IV. 18) 

yfl - e 2 sin 2 (Ai) 

Consider a Cartesian coordinate system with its origin at the center of the 
WGS84 ellipsoid, oriented such that its z-axis is aligned due north, its x-axis intersects 
the prime meridian, and its y-axis completes the right hand rule. This reference 
coordinate system is termed Earth Centered Earth Fixed (ECEF). Then, given the 
position reported by GPS (IV.17), the position expressed in the ECEF reference frame 
is 



x = (N + h) cos(A 2 ) cos(Ai). 
y = (jV + h) cos(A 2 ) sin(Aj ). 

2 = (iV(l - e 2 ) + h) sin(A 2 ). 

Let [Aj 0 A 2o h] T be the surveyed position of the workstation, and let [xo, yo, 2o] T 
the same location expressed in ECEF coordinates, and suppose this position is used 
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to define the origin of a local tangent plane reference frame. Sign convention and 
orientation of the local tangent plane reference frame are defined as positive x-axis 
values extending due north, positive y-axis values extending due east, and positive 
z-axis values extending straight down. Furthermore, assume [xf Tog yf TO g z frog] T is the 
reported position of the vehicle expressed in the ECEF coordinate system. Then, 



Xl 




-sin(A 2o )cos(A, 0 ) - sin(A 2o ) sin(A lo ) cos(A 2o ) 




(Zfrog %o) 


yi 


— 


— sin(Aj 0 ) cos(A lo ) 0 




( Hfrog ~ Vo) 


Zi 




cos(A 2q ) cos(A] 0 ) -cos(A 2o )sin(Ai 0 ) -sin(A 2o ) 




( z Jrog z o) 



(IV. 19) 



is the position of the vehicle in the local tangent plane reference frame. This position 
was displayed on the console, and later, the trajectories tracked by Frog were defined 
in this reference frame. 

The final issue that needed to be resolved was the accurate and timely calibra- 
tion of the command signals from the console to the vehicle. The autopilot onboard 
the vehicle expected to see commands in terms of PWM signals. The pulse width 
modulated signals were generated by first converting the command signals to analog 
voltages via the Digital-to-Analog converter. Then the analog voltages were directed 
to a Futaba transmitter, specially modified to respond to analog voltages instead of 
the movement of the exterior joysticks. Next, a receive module onboard Frog con- 
verted the transmitted signals to PWM format, where they were sent to the autopilot. 
A functional block diagram of the architecture used to calibrate the uplink is shown 
in Figure 52. 

The key to calibrating the uplink was the determination of the functions F~J m 
(Block l) and G~l lt (Block 2) in Figure 52. This was done in two steps. The first 
step involved sending constant commands to the DAC from the user interface (Block 
10, 11 and 12). This resulted in flight patterns consisting of numerous constant 
rate turns to the left and to the right, and to constant rate climbs and descents. 
A receive module, which was a duplicate of the one onboard Frog , converted the 
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Figure 52. Calibrating the Uplink 



transmitted signal to PWM format. This was then converted and stored in digital 
format by the RFTPS via the pulse width modulated signal converter I/O module 
(Block 3, 4 and 5). Concurrently, the onboard navigation sensors transmitted yaw 
rate and climb rate data, among other things, to the RFTPS (Block 6, 7, 8 and 9). 
Post processing of the data revealed linear relationships between the PWM signals 
sent to the autopilot and the respective climb/descent or turn rate of the vehicle 
(see Figure 53). The relationships between the command signals in terms of their 
PWM values (r Cpwm and h Cpwm ) and the steady-state performance of Frog in flight are 
expressed as 



0.0885r Cpu , m - 151.476, 


(IV. 20) 


piD m ( ^Cpxn m ) ? 




12.3929A Cpwm -18193, 


(IV. 21) 


^ 'h pW m (^ c pu/m)? 


(IV. 22) 
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where r is the vehicle’s yaw rate in degrees per second, and h is the vehicle’s climb/descent 
rate in feet per minute. 





PWM 

Figure 53. Top graph: measured vehicle yaw rate (degrees per second) versus PWM 
signal (// seconds) into the lateral channel of the inner-loop autopilot. Bottom graph: 
measured vehicle climb rate (feet per minute) versus PWM signal (//seconds) into the 
longitudinal channel of the inner-loop autopilot. 

The second step in calibrating the uplink involved the modified Futaba trans- 
mitter (Block 2). The known voltages applied to the transmitter were recorded, 
while the duplicate receive module was used to decode the transmitted signal, and 
the RFTPS was used to record the corresponding PWM value (Block 3, 4 and 5). 
The relationships were linear and expressed as 



r Cpwm = 686.7r Wt - 234.1, 


(IV. 23) 


= - ^T v0l X r V0lt)l 




h C P u,m = 800 h vo it - 470, 


(IV.24) 
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— : G kvolt (Kou), 

where h vo n is voltage sent to the longitudinal channel of the modified Futaba trans- 
mitter, and r vo it is the voltage sent to the lateral channel. 

That completed the calibration of the uplink. Given a commanded climb/descent 
rate ( h c ) in feet per minute, and commanded yaw rate (r c ) in degrees per second at 
the console, the voltages sent to the modified transmitter were calculated as 

Kelt = Gf Fr 1 h c , (1V.25) 

n volt H'p wrn v ' 

Tvoit = G T * o[t F rp l wm r c . (IV. 26) 

E. RESULTS AND ANALYSIS 

This section summarizes the results of two flight tests of the guidance and 
control algorithm developed in the preceding section. While the trajectory definition 
parameter was identical for both tests, the test flights took place on different days 
with very different environmental conditions. Furthermore, on the first test, only the 
lateral channel of the guidance algorithm was active. Longitudinal control was done 
open-loop by the console operator. Both flights took place at approximately 500 feet 
of altitude above ground level. Wind measurements were made at ground level on 
the runway, and later incorporated into the simulation. The pilot maintained control 
of the throttle throughout the tests. During autonomous flight, the pilot left the 
throttle at the trim setting for the trajectory defined. The rudder was not used. 

At the time of the tests, the RFTPS had been in operation for approximately 
three months. The two flight tests presented were not the first time the RFTPS 
was used to control Frog. Some of the model identification process, and most of 
the calibration process, used the RFTPS to control Frog. This was primarily due 
to the fact that the pilot used a conventional stick configuration to control Frog. 
Longitudinal signals corresponded to fore and aft motion, and lateral commands 
corresponding to left and right motion. Movement along one axis invariably corrupted 
the signal along the other axis. In contrast, control signals generated by the RFTPS 
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were very precise and perfectly decoupled. Also, tests had been successfully completed 
whereby a portable computer was used to send commands to the RFTPS and control 
Frog. The flight tests presented do represent the first time that a guidance and control 
algorithm was used to control Frog fully autonomously. 

1. First Flight Test 

The first test of the guidance algorithm occurred on June 20, 1997. The 
wind was measured at 10 mph out of the northwest with gusts to 15 mph. The 
pilot launched Frog and brought it up to a safe altitude. Control was passed to the 
console operator who turned the controller on. The trajectory definition is repeated 
for convenience, 



7 C = 0 degrees, 

i[) c = 5 degrees per second, 

b c — 1146 feet. 

Autonomous control of the longitudinal channel was not implemented. The GPS 
operated in the differential position mode throughout the test. Prior testing of this 
GPS revealed that position error in the horizontal plane, using GPS in differential 
mode, had a standard deviation of 10 feet. The positions shown are those reported 
by GPS. The true position is unknown but thought to be within about 10 feet of that 
shown. Later, the process was repeated in simulation using the measured value of 
the wind. The trajectory flown by the vehicle is shown projected onto the horizontal 
plane in Figure 54. Additionally, the reference trajectory and the trajectory obtained 
in simulation is shown. 

The controller activity in the lateral channel is shown alongside the lateral 
error signal in Figure 55. According to the orientation of the Frenet frame, the error 
signal went negative when the vehicle was inside of the reference trajectory. The 
controller activity has noticeable high frequency content caused by the interaction 
of the one second update rate of the GPS and the slow complex zeros placed in 
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Figure 54. The first (near) autonomous flight of Frog. The commanded trajectory is 
shown as a solid line while the path flown is shown as a broken line. Also shown are the 
results from nonlinear simulation. The wind was 18 feet per second out of the west, 
northwest. The graph is oriented with the y-axis pointing due north. Approximately 
two and one half minutes of flight are shown. 



the controller. The nominal commanded yaw rate should be 5 degrees per second, 
but instead has a mean value of 8 degrees per second. This is indicative of the 
integral action compensating for modeling and calibration errors. The ground speed 
(not shown) oscillated by twice the magnitude of the measured wind. The vehicle’s 
indicated airspeed was nearly constant. Therefore, the wind appeared as a sinusoidal 
disturbance at a frequency of 0.796 radians per second. As the loop gain cross over 
frequency was 0.5 radians per second along the lateral channel, little gain was available 
to suppress this disturbance. 

Similar qualitative and quantitative performance was observed in simulation. 
The slight difference in orientation of the two trajectories is most likely due to mis- 
alignment of the modeled and actual wind. Table VI quantitatively compares flight 
test data to linear and nonlinear simulation results. The close agreement in average 
(RMS) tracking error indicates that the design has good robustness properties. The 
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Figure 55. Top graph: Commanded vehicle yaw rate in degrees per second from the 
controller along the trajectory, flown on June 20, 1997. Bottom graph: Lateral error 
in feet, as computed by the guidance algorithm, along the same trajectory. 



relatively large tracking errors are indicative of the low loop gain available for dis- 
turbance rejection of wind effects. In retrospect, performance requirements should 
be formulated to ensure adequate gain at these frequencies. The frequency of this 
disturbance, u>d, is well known at the design stage. It can be defined using fj c as 






360 

i>c cos ( 7 C ) 



(IV. 27) 



With estimates of the maximum magnitude of the wind, and requirements on the 
maximum permissible tracking errors allowed, a performance requirement could nat- 
urally be formulated as an constraint. 
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Linear Simulation 


Nonlinear Simulation 


Flight Test 


Lateral Mean Error 


-21.92 feet 


-22.1 feet 


-18.35 feet 


Lateral RMS Error 


111.27 feet 


111.49 feet 


121.24 feet 


Lateral Peak Error 


140.57 feet 


140.01 feet 


219.54 feet 



Table VI. Error Comparison for Flight of June 20, 1997: Flight Test & Simulation 



2, Second Flight Test 

On July 23, 1997, a second flight test was conducted. The second flight test 
provided the opportunity to test the guidance algorithm in near ideal weather condi- 
tions. The wind was 1 to 2 miles per hour out of the west, and the air was absent of 
thermal activity and disturbances. This time, the longitudinal channel was operative. 

Frog was brought up to altitude by the pilot and control handed over to the 
console operator. The console operator switched to autonomous flight with the tra- 
jectory parameters defined as on June 20, 1997. Once again, post test flight data is 
compared with simulation results incorporating wind levels measured. The results 
are summarized in Table VII, followed by data from the flight. 





Linear Simulation 


Nonlinear Simulation 


Flight Test 


Lateral Mean Error 


-2.419 feet 


-2.431 feet 


5.49feet 


Lateral RMS Error 


11.162feet 


11.1869 feet 


21.41 feet 


Lateral Peak Error 


15.24 feet 


15.21 feet 


41.5 feet 


Longitudinal Mean Error 


-.006 feet 


-.005 feet 


3.658 feet 


Longitudinal RMS Error 


.058 feet 


.091 feet 


14.949 feet 


Longitudinal Peak Error 


.089 feet 


.0584 feet 


37.1725 feet 



Table VII. Error Comparison for Flight of July 23, 1997: Flight Test h Simulation 



The resulting path in space flown on July 23 is shown in Figure 57. It is 
compared to the reference trajectory. The projection of the trajectory tracked onto 
the horizontal plane is shown in Figure 56. It, also, is compared to results from 
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simulation. The history of the controller activity is shown for the lateral channel in 
Figure 58, and for the longitudinal channel in Figure 59. 




Figure 56. The second autonomous flight of Frog. Weather conditions were near 
ideal as Frog tracked the 3-D trajectory shown by the dashed line. The results from 
simulation are shown as a dash-dot line. 




Figure 57. The second autonomous flight of Frog. Weather conditions were near ideal 
as Frog tracked the 3-D trajectory shown. 
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Figure 58. Top graph: Commanded vehicle yaw rate, in degrees per second from the 
controller, along the trajectory flown on July 23, 1997. Bottom graph: Lateral error 
in feet, as computed by the guidance algorithm, along the same trajectory. 



With little wind, the test data characterizes the performance of the guidance 
and control algorithm in the presence of modeling errors, transport lag, hardware non- 
linearities, and sensor noise. The simulation results were a close match to flight test 
data along the lateral channel. The average tracking error just over 21 feet. It can be 
seen in Figure 58 that the integral control is holding about one half degree per second 
commanded yaw rate to compensate for constant disturbances. In flight, tracking 
errors along the longitudinal channel were even less than lateral tracking errors. The 
average command signal along the longitudinal channel was —5 feet per second. This 
probably indicates that the power setting was too high. Subsequent testing will 
incorporate airspeed and throttle control which will address this issue. In simulation, 
however, longitudinal tracking errors were extremely small. The discrepancy is most 
likely the result of unmodeled vertical disturbances in the airmass due to light thermal 
activity. 
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Figure 59. The output of the longitudinal channel of the controller is shown along 
with the error signal on which it is acting. Note that the positive z-axis points down 
and that positive error corresponds to the vehicle being below the reference trajectory. 

F. CONCLUSIONS 

In this chapter, an integrated system for the design, development and testing of 
guidance, navigation, and control algorithms for unmanned air vehicles was presented. 
Extensive use was made of commercial off-the-shelf (COTS) hardware. This kept costs 
low, and made the system easily scaleable. Sophisticated autocode tools allowed a 
two man team to write, test, and maintain thousands of lines of error free real time 
code. In a single, unified environment, the avionics system was designed, simulated, 
simulated incorporating hardware-in-the-loop (HITL), tested in flight, and used to 
control an aircraft in flight. 

As a proof-of-concept demonstration, the RFTPS was used to take the in- 
tegrated guidance and control concept presented in Chapter III and evaluate it in 
the environment for which it was proposed to he used. The project began with an 
unmanned air vehicle with unknown flight characteristics. In addition, a “black box” 
autopilot with unknown internal dynamics was placed onboard the vehicle. Through 
the use of the RFTPS, a high fidelity simulation of the vehicle and autopilot was built 
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and used to synthesize the applicable guidance and control laws. Extensive testing 
in simulation followed, and appropriate user interfaces were developed. Then, the 
RFTPS was used to control the vehicle in flight using the guidance and control laws 
developed, collect the data during the flight test, post process the data, and evaluate 
the performance of the algorithms. 

The success of the project demonstrated both the utility of the integrated guid- 
ance and control algorithm as well as the capabilities of the RFTPS. The integrated 
guidance and control algorithm was shown to work well in a real world application. 
Performance in flight was very close to performance in simulation, which speaks well 
of the robustness properties of the controller. Changing environmental conditions 
highlighted some disturbance rejection issues that should be addressed in the future. 
The RFTPS was shown to be powerful, portable, rugged, effective in the field and in 
the lab, and safe and reliable at controlling an aircraft. 



137 







138 




APPENDIX A. THE TAIL SIZING DESIGN 

TOOL 



1. INTRODUCTION 

This appendix details the operation of the Tail Sizing Design Tool The Tail 
Sizing Design Tool runs from the MATLAB command line, and requires the LMI and 
CONTROL Toolboxes. The purpose of the Tail Sizing Design Tool is to allow the 
user to map and view the Tail Sizing Design Space for a given aircraft definition, or 
compare the Tail Sizing Design Space of two different aircraft definitions. There are 
three versions of the Tail Sizing Design Tool The Tail Sizing Design Tool Version 
A maps the Tail Sizing Design Space for the high angle-of-attack excursion search- 
ing over static, state feedback, controllers. Version B solves the same problem but 
searches over dynamic, output feedback, controllers. Version C maps the Tail Sizing 
Design Space for the wind-shear penetration constraint searching over static state 
feedback controllers. Version A and B are well developed and allow the user to in- 
clude aeroelastic effects, the addition of a canard, and to view the Tail Sizing Design 
Space in a three dimensional plot. Version C is less well developed and provides for 
only rigid body, tail only, aircraft dynamics, and two dimensional views. 

The general structure of all three versions is similar. The Tail Sizing Design 
Tool is comprised of three directories which must reside on the same level within 
the user’s file system. The first directory is termed (7W, which is an acaronym for 
graphical user interface. As the name suggests, it contains files which allow the 
user to specify an aircraft definition, generate a database for an aircraft definition or 
view a Tail Sizing Design Space for an aircraft definition by pushing the appropriate 
“buttons” on the display generated. The second directory is termed MODEL , and 
it contains the files necessary to model the dynamics of an aircraft. This directory 
contains the user supplied aerodynamic stability and control derivatives. The third 
directory is termed DATA , and it contains data files generated by running the Tail 



139 



Sizing Design Tool. 

Assume for the moment that the proper stability and control derivatives have 
been placed appropriately in the MODEL directory. The process is as follows. 

1. From the GUI directory, run size.it. 

2. Fill in the appropriate values for the range of tail volumes and actuator rates 
of interest. 

3. Specify the aircraft definition parameters such as the use of a canard or inclu- 
sion of a flexible mode. 

4. For output feedback controllers, specify the output sensors to be used from 
the available choices. 

5. Name the database file and start the generation of the data. This executes 
the program getdata. 

6. Run the program view.it to querry a database previously generated. 



2. GENERATING A DATABASE 

First, we will consider the two programs size.it and getdata. The interface 
generated by size.it (Version B) is shown in Figure 60. By selecting the push button, 
” Output Page,” another GUI is created (Figure 61). This provides the user the option 
of specifying which sensors to use for feedback. When the user is satisfied with the 
definition of the aircraft, he begins execution by selecting the pushbutton Run It!. 
This executes the program Getdata , shown next for Version B. 



raramramramrarara 



Getdata. m 



mmmmmrammmra 



7 , Initialize global variables and move to the model directory. 
cgrate=[] ; data=[]; 
cd . . 
cd model 

‘/.Use the graphics handles to transfer the user defined aircraft 
*/,def intions to variables in the workspace. Begin with the frequency 
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Figure 60 . Interface called by size Jt. in 



'/, of the flexible mode. 

freql = str2num(get (freq_ed, ’ string' )) ; 
save flexdata freql 

'/.This is the location of the flexible mode sensors. 
sens_loc = str2num(get(sens_ed, 'string')) ; 

'/.Some locations (nodes/anti-nodes) make the problem 

'/.unobservable. Check for this and alert the user if it is a problem, 
[phi , phidot] =modeshpe (sens_loc'/,250/ 100) ; 
if (phidot < le-5) , 

disp('Flexible pitch rate is nearly unobservable.') 
disp('Change the sensor location ... stopping’) 
elf 

'/, Message 

frame_title = uicontrol (gcf , . . . 

'style' , 'text' , . . . 

'position' , [120 270 200 30] .... 

'string' , 'There was a problem, ’ , . . . 

’foregroundcolor [1 1 1],... 

'backgroundcolor ’ , [0 00]); 

'/, Message 

frame_title = uicontrol (gcf , . . . 
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Figure 61 . Selection of output sensors 



’ style ’ , ’ text ’ , . . . 

'position’ , [120 240 200 30] , . . . 

'string' , 'see the command window. ' , . . . 

’foregroundcolor ’, [1 1 1],... 

’ backgroundcolor’ , [0 00]); 
clear 
break; 
end; 

‘/.Name the file to store the data, 
opfname = get(save_ed, ’string’) ; 

'/.Integer 1 if model is aeroelastic. 
flexflag = get (ck_f lex , ’value’ ) ; 

'/.Integer 1 if the model has a canard, 
canardflag = get (ck_canard, ’value’) ; 

'/.If a canard is present, specify its canard volume, 
if (canardflag==l) , 

split = str2num(get (canard_ed, ’ string’ )) ; 
else, split = 0; end; 

'/.This is the range and increments of tail volume to check, 
vhlo = str2num(get (minvol_ed, ’ string’ )) ; 
vhhi = str2num (get (maxvol.ed, ’ string’ )) ; 
vhinc = str2num(get (incvol_ed, ’ string’ )) ; 

‘/.This is the required closed loop damping, 
damping = str2num(get (damp.ed, ’ string’ )) ; 
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theta = 2'/,acos (damping) ; */ pole placement region definition 
'/.This is the range and increment of actuator rates, 
ratelo = str2num(get (minrate_ed , ' string' )) ; 
ratehi = str2num(get(maxrate_ed, 'string')) ; 
rateinc = str2num(get(incrate_ed, ’ string' )) ; 

'/This specifies which sensors to use for output feedback. 

'/ An integer 1 indicates use of the sensor. Choices are 

'/. flight path angle, airspeed angle of attack, local 

'/ pitch rate, local pitch angle. The default is full 

'/, information output feedback. 

gamaflag = get (ck_gama, 'value' ) ; 

asflag = get (ck_as, 'value') ; 

aoaflag = get (ck_aoa, ' value ') ; 

lprflag = get (ck_lpr, 'value' ) ; 

lpaflag = get (ck_lpa, 'value ') ; 

fsflag = get (ck_fs, 'value') ; 

% The actuator dynamics are specified in the function 
% ACTAUTOR.M. Canard and tail actuators do not have 
7 , to be the same although dummy dynamics must be 
•/ specified for a canard even if not used. 

[Aa , Ba , Ca , Da] = actuator; 
act=length(Aa) ; 

'/. For the first order actuator modelused here, 

'/. specify the time constant 
a_bw = 20; 

'/, The angle of attack excursion is equal to acos(gust/vt) . 
gust = 66; /fps 

'/, The equilibrium point is specified by the flight 
7 . path angle and true airspeed of the HSCT. 
gamma = 0; vt = 422; '/ deg fps 

'/. Specify tolerance limit for the eg binary search, 
tol = 0.01; 

'/, MATLAB uses the matrix variable ''region'' to define 
'/ the closed loop poleplacement region. Form it 
/ based on user inputs . 
region=zeros (4 , 8) ; 
region (1, l)=.02+li; 
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region (2, 2)=-40+li ; 
region(3,3)=0+2i ; 
region(l,5)=l; 
region(2,6)=-l; 
region(3,7)=cos(damping) ; 
region(4,8)=cos(damping) ; 
region(3 , 8) =-sin (damping) ; 
region(4,7)=sin(damping) ; 

*/, Set up three loops to map the design space. 

'/, First, loop through range of tail volumes, 
for vh=vhlo :vhinc : vhhi ; 

% Initialize some variables for each tail volume. 
yemax=0; ydemax=0; ycmax=0; ydcmax=0; 

*/, Loop through the range of actuator rates, 
for ratelimit = ratelo :rateinc:ratehi; 

*/, Reset eg upper and lower bounds for the 

'/. binary search. 

xegmax = 1.0; xegmin = 0 ; 



% Binary search on c.g. ;drop out when 

y, within tolerance. 

while (norm(xcgmax-xcgmin) > tol) ; 

save cgvh xcg vh split; 

'/, Load appropriate data for trim; initial 
y, guesses for trim, 
load trimdata 

*/, Trim at eg specified and linearize at peak 
'/, of gust condition, 
if ( flexflag == 1) , 

[xtrim,utrim,ytrim] = trim( , eom_b70w' ,x0,u0, . . . . 
[vt ; gamma] , [] , [1] , [1 ; 2] ) ; 
xtrim(2) = xtrim(2) + gust; 

[af ,bf , cf ,df] = linmod( ’eom.byOw’ ,xtrim,utrim) ; 
bf = bf (: ,1:2) ; df = df (: ,1:2) ; 
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xtrim(2) = xtrim(2) - gust; 
else , 

[xtrim,utrim,ytrim] = trim( ’eom_b70r’ ,xOr,uO, . . . . 

[vt ; gamma] , [] , [l] , [1 ;2] ) ; 
xtrim(2) = xtrim(2) + gust; 

[af ,bf ,cf ,df] = linmod( ' eom_b70r' , xtrim,utrim) ; 
bf = bf (: ,1 :2) ; df = df(:,l:2); 
xtrim(2) = xtrim(2) - gust; 
end; 

*/, Connect the actuators in series. 

[A,B,C,D] = series(Aa,Ba,Ca,Da,af ,bf ,cf ,df) ; 
[m,n]=size(A) ; 

7 . Place the actuator states after the plant states, 
[n ,p] =size(B) ; 

T1 = [zeros(n-2,2) eye(n-2) ; eye(2) zeros(2,n-2)] ; 

A = TT/.A/T1; 

B = T1*/.B; 

C = C/Tl ; 

*/, Similarity transformation used to make local 
*/, pitch rate and local body angle states vice 
'/. using generalized coordinates, 
if (flexf lag==l) , 

[phi , phidot] =modeshpe(sens_loc'/,250/ 100) ; 

T2 = [eye (4) zeros (4, 2+act) ; 

[0001 phidot 0] zeros (1, act) ; 

[00100 phidot] zeros ( 1 , act) ; 
zeros(act,6) eye(act)] ; 
else, T2=eye(n); end; 

Abar = T2A/T2; 

Bbar = T2B; 

Cbar = C/T2 ; 

'/, Form output matrix based on user specified sensors 
'/, airspeed is Cbar(l,:) 

'/, flight path angle is Cbar(2,:) 

*/. angle of attack is Cbar(3,:) 

*/, local pitch rate is 5th state [0 0 0 0 1 0 0 0] 

# /, local pitch angle is 6th state [0 0 0 0 0 1 0 0] 
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C_select= [] ; 
if gamaf lag == 1 ; 

C.select = [C_select ;Cbar(2 , : )] ; 
end; 

if asflag == 1; 

C_select = [C_select ;Cbar(l , : )] ; 
end; 

if aoaf lag == 1 ; 

C_select = [C_select ;Cbar(3 , : )] ; 
end; 

if lprflag == 1; 

C_select = [C_select ; [zeros (1 ,n-4) 100 0]]; 
end; 

if lpaf lag == 1 ; 

C_select = [C_select ; [zeros (1 ,n-4) 010 0]]; 
end; 

if fsflag == 1; 

C_select = eye(n); 
end; 

Cbar=C_select ; 

[n,m] =size(Bbar) ; 

[p ,n] =size(Cbar) ; 

'/, Use the angle of attack excursion to define 
‘/, the initial condition. 

i0=[0 gust 0 0 zeros(l,n-4-act) zeros (1, act)] ’ ; 
*/,i0=T'/,i0 ; 

*/, Specify rate and ampltude constraints. The 
'/, rate limit comes from the 

*/, GUI. The saturation limit is set at 30 degrees. 



ucmax = 30 ; uemax = 30 ; 

udcmax = ratelimit; udemax = ratelimit 



'/, Finally, form the LMIs. 
setlmis ( [] ) ; 

R = lmivar(l, [n l] ) ; '/.R*/, 

S = lmivar(l , [n l] ) ; ‘/. S '/. 
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Ak = lmivar(2, [n n] ) ; '/, Ak '/, 

Bk = lmivar(2, [n p] ) ; % Bk ’/, 

Ck = lmivar(2, [m n] ) ; "/, Ck */, 

’/, This is pole placement LMIs formed using 
’/, LMI region matrix REGION. 

[rr,rc]=size(region) ; 
if rc~=2‘/rr, 

error ( 'REGION must be a Mx(2M) matrix'); 

end 

rl=0 ; 

’/, scan each region (L,M) (diagonal blocks) 
while rl<rr, 

bs=round(imag(region(rl+l ,rl+l) ) ) ; 
if ~bs, bs=rr-rl; end 

L=real(region(rl+l :rl+bs,rl+l :rl+bs) ) ; 
M=region(rl+l : rl+bs,rr+rl+l : rr+rl+bs) ; 
pole=newlmi; “/, pole placement in the region (L,M) 
for i=l:bs, 

nbi=2'/,i-l; '/, row of block in target LMI 

'/, off-diagonal 2x2 block 
for j=l : i-1 , 

nbj=2'/,j-l; '/, col of block in target LMI 
lmiterm( [pole ,nbi ,nbj ,R],L(i,j),l); 
lmiterm( [pole ,nbi ,nbj ,R] , Abar,M(i , j ) ) ; 
lmiterm( [pole ,nbi ,nbj ,R] , M ( j , i) , Abar ' ) ; 
lmiterm( [pole ,nbi ,nbj ,Ck] ,Bbar,M(i, j)) ; 
lmiterm( [pole ,nbi ,nbj , -Ck] ,M(j , i) ,Bbar ' ) ; 
lmiterm( [pole,nbi+l ,nbj+l ,S] ,L(i , j) , 1) ; 
lmiterm( [pole ,nbi+l ,nbj+l ,S] ,M(i , j ) , Abar) ; 
lmiterm( [pole ,nbi+l ,nbj+l ,S] ,Abar ' ,M(j , i) ) ; 
lmiterm( [pole ,nbi+l ,nbj+l ,Bk] ,M(i , j ) ,Cbar) ; 
lmiterm( [pole ,nbi+l ,nbj+l ,-Bk] ,Cbar ' ,M(j , i) ) ; 
lmiterm( [pole ,nbi+l ,nbj , 0] ,L(i , j) ) ; 
lmiterm([pole,nbi+l,nbj ,Ak] ,M(i, j) , 1) ; 
lmiterm([pole,nbi+l,nbj ,0] ,M(j ,i)7.Abar') ; 
lmiterm( [pole ,nbi ,nbj+l ,0] , L(i , j ) ) ; 
lmiterm( [pole ,nbi ,nbj+l , -Ak] ,M(j , i) , 1) ; 
lmiterm( [pole ,nbi ,nbj + l ,0] ,M(i , j )*/,Abar) ; 
end 

V. diagonal 2x2 block 

lmiterm( [pole,nbi,nbi ,R] ,L(i, i) , 1) ; 
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ImitermC [pole ,nbi ,nbi ,R] ,Abar ,M(i, i) , ' s' ) ; 
ImitermC [pole ,nbi ,nbi ,Ck] ,Bbar ,M(i , i) , ’ s ' ) ; 
ImitermC [pole,nbi+l ,nbi+l ,S] ,L(i ,i) , 1) ; 
lmiterm( [pole,nbi+l ,nbi+l ,S] ,M(i ,i) ,Abar , ' s ’ ) ; 
ImitermC [pole,nbi+l ,nbi+l ,Bk] ,M(i , i) ,Cbar, ’s'); 
ImitermC [pole,nbi+l ,nbi ,0] ,L(i ,i)) ; 
lmiterm( [pole,nbi+l ,nbi ,Ak] ,M(i ,i) , 1) ; 

ImitermC [pole ,nbi+l ,nbi , 0] ,MCi , i)'/,Abar ’ ) ; 
end 

rl=rl+bs ; 
end ’/ while 



’/, This is x_0 LMI from the invariat elipsoid. 

x0_lmi = newlmi ; 

lmiterm C [~x0_lmi 1 1 0], 1); 

ImitermC [~x0_lmi 210], iO); 

ImitermC [-x0_lmi 2 2 R],l,l); 

ImitermC [-x0_lmi 3 1 S] , 1 , iO) ; 

ImitermC [-x0_lmi 320], eyeCn)); 

ImitermC [-x0_lmi 3 3 S] , 1, 1); 

'/, This is Cu_max) saturation limit elevator LMI. 
umaxe_lmi = newlmi ; 

ImitermC [~umaxe_lmi 1 1 R] , 1,1); 

ImitermC [-umaxe.lmi 2 1 0], eyeCn)); 

ImitermC [~umaxe_lmi 2 2 S] , 1,1); 

ImitermC [~umaxe_lmi 3 1 Ck] , [0 1] , 1) ; 

ImitermC [~umaxe_lmi 3 2 0], 0'/,eyeCn)); 

ImitermC [~umaxe_lmi 330], uemax~2) ; 

•/. This is (u_max) saturation limit canard LMI. 
if Ccanardf lag==l ) , 
umaxc_lmi = newlmi; 

ImitermC [-umaxc_lmi 1 1 R] , 1,1); 

ImitermC [-umaxc_lmi 2 1 0],eye(n)); 

ImitermC [~umaxc_lmi 2 2 S] , 1,1); 

ImitermC [-umaxc_lmi 3 1 Ck] , [l 0],l); 

ImitermC [~umaxc_lmi 3 2 0], O'/.eyeCn)); 

ImitermC [~umaxc_lmi 330], ucmax~2); 
end; 
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'/, This is (ud_max) actuator rate elevator LMI . 

udmaxe.lmi = newlmi; 

lmiterm( [-udmaxe_lmi 1 1 R] , 1,1); 

lmiterm( [-udmaxe_lmi 2 1 0],eye(n)); 

lmiterm( [~udmaxe_lmi 2 2 S] , 1,1); 

lmiterm( [-udmaxe_lmi 3 1 Ck] , [0 a_bw],l); 

lmiterm( [-udmaxe_lmi 3 1 R] , [zeros(l ,n-l) -a_bw],l); 

lmiterm( [-udmaxe_lmi 320], [zeros(l ,n-l) -a_bw] ) ; 

lmiterm( [-udmaxe_lmi 330], udemax~2) ; 

X This is (ud_max) actuator rate canard LMI . 

if (canardflag==l) , 

udmaxc_lmi = newlmi ; 

lmiterm( [-udmaxc_lmi 1 1 R] , 1,1); 

lmiterm( [-udmaxc.lmi 2 1 0],eye(n)); 

lmiterm( [-udmaxc_lmi 2 2 S] , 1,1); 

lmiterm( [-udmaxc_lmi 3 1 Ck] , [a_bw 0],1); 

lmiterm( [-udmaxc_lmi 3 1 R] , [zeros(l ,n-2) -a_bw 0],l); 

lmiterm( [~udmaxc_lmi 320], [zeros (l ,n-2) -a_bw 0]); 

lmiterm( [-udmaxc_lmi 3 3 0], udcmax~2) ; 

end; 

lmisys = getlmis; 

'/, Chech to see if the set of LMIs is feasible. 

X TMIN should be negative for feasibility and 
*/, XFEAS is the decision vector that validates 
*/, the above LMIs. 

[tmin,xf eas] = feasp(lmisys, [0 200 le8 0 0] ) ; 

’/, Values of TMIN less than 0 are not strictly feasible 
*/, but seem to be close enough to work, 
if tmin < le-4; 

*/, Recover the LMI variables. 

R = dec2mat (lmisys, xf eas, R) ; 

S = dec2mat (lmisys, xf eas, S) ; 
cik = dec2mat (lmisys ,xf eas ,Ak) ; 
bk = dec2mat (lmisys ,xf eas ,Bk) ; 
ck = dec2mat (lmisys ,xf eas ,Ck) ; 
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'/, Construct the controller and form the close loop system 
[u,sd,v]=svd(eye(n)- R'/,S) ; */, factorize I-RS 

sd=diag(sqrt (1 ./diag(sd))) ; 

Ni=sd*v'; Mti=u*sd; 

Ac=Ni*(ak-S* (Abar)R-bkCbarR-S*Bbar*ck)*Mti ; 

Bc=Ni* (bk) ; 

Cc=(ck) *Mti ; 

P = ltisys (Abar , Bbar ,Cbar,zeros(p,m)) ; 
cont = ltisys (Ac ,Bc ,Cc,zeros(m,p) ) ; 
clsys = slf t (P , cont ,m,p) ; 

[acl,bcl,ccl,dcl] = ltiss(clsys) ; 

'/, Simulate the closed loop system with the specified 
*/ initial condition and record the maximum actuator 
*/ 0 amplitudes and rates. 

[ye,x,t] = initial (acl, [Bbar(: ,2) ;zeros(n,l)] ,... . 
[zeros(l,n) [0 1] *Cc] ,0, [i0;zeros (n, 1)] ) ; 
yemax = max(abs (ye) ) ; 

[yc,x,t] = initial(acl , [Bbar( : , 2) ;zeros(n, 1)] , . . . . 
[zeros(l,n) [1 0] *Cc] , 0 , [iO ; zeros (n, 1)] ) ; 
yemax = max(abs(yc)) ; 

[yde,x,t] = initial (acl , [Bbar( : ,2) ;zeros(n, 1)] . 

[ [zeros (1 ,n-l) -a_bw] [0 a_bw]*Cc] ,0, [i0;zeros(n,l)]) ; 
ydemax = max(abs(yde)) ; 

[ydc,x,t] = initial (acl , [Bbar( : ,2) ;zeros(n, 1)] . 

[ [zeros (1 ,n-2) -a_bw 0] [a_bw 0] *Cc] , 0 , [iO ;zeros (n , 1)] ) ; 
ydemax = max(abs(ydc)) ; 
eigvalue = eig(acl) 

data = [data; xcg ratelimit vh split yemax yemax ydemax.. 
ydemax tmin utrim' xtrim’ eigvalue']; 

'/ Move eg limit aft. 
xegmin = xcg; 

X Otherwise, move eg limit forward, 
else, 

xegmax = xcg; 

end; ’/SET FEASIBLE, if 

*/, Chech to see if aft eg limit is 
’/, within tolerance 

if ( abs (xegmax-xegmin) < tol | xcg > .95 | xcg < .05 ) 
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break ; 

end; °/ 0 WITHIN TOLERANCE, if 
end; ‘/.BINARY CG SEARCH, while 

°/. This is the aft most eg for the rate, tail 
% volume combination. 
if(ydcmax ~= 0 I ydemax ~= 0), 

cgrate= [egrate ; vh xcg ydemax yemax ydemax yemax eigvalue']; 
end; 

end; ‘/.RATE SWEEP, for 
end; ‘/.VH SWEEP, for 

'/. Change to the DATA directory and store the 
‘/ P result with the name provided from the GUI. 

'/, Return the user to the GUI directory, 
cd . . 
cd data 

eval(['save * , opf name, > .mat egrate data ’ ] ) 
cd . . 
cd gui 



The above code is associated with the output feedback problem. For the state 
feedback problem, the code is very similar. The essential difference is in the LMI 
formulation and reconstructionof the controller. That part of the code from Version 
A substantially different from Version B is shown next. 

Form the state feedback LMIs 

setlmis ( [] ) ; 

Y = lmivar(l,[n 1]); 

W = lmivar(2,[2 n] ) ; 

*/, This is Lyapunov lmi with real part constraint 
lyap = newlmi ; 

lmiterm( [lyap 1 1 Y] , Abar+center*eye(n) , 1, ; s ; ); 
lmiterm( [lyap 1 1 W] , Bbar, 1, ; s ; ); 

This is the pole placement LMI 
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pplac = newlmi ; 

lmitermC [pplac 1 1 Y] , sin(theta/2)*Abar, 1, ’s’); 

lmitermC [pplac 1 1 W] , sin(theta/2)*Bbar, M, 's’); 

lmiterm( [pplac 2 1 Y] , cos (theta/2) *Abar, 1); 

lmitermC [pplac 2 1 W] , cos(theta/2)*Bbar, M) ; 

lmitermC [pplac 2 1 -Y] , -cos (theta/2) , Abar’) ; 
lraiterm( [pplac 2 1 -W] , -cos(theta/2)*M’ ,Bbar’ ) ; 
lmitermC [pplac 2 2 Y] , sin(theta/2)*Abar , 1, ’s’); 

lmiterm( [pplac 2 2 W] , sin(theta/2)*Bbar , M, ’s’); 

V. This is the Y > 0 LMI 

ypos = newlmi; 

lmitermC [-ypos 1 1 Y],l,l); 

*/„ This is x_0 LMI 

xO_lmi = newlmi ; 
lmitermC [~xO_lmi 1 1 0] , 1) ; 
lmitermC [-x0_lmi 2 10], iO) ; 
lmitermC [-x0_lmi 2 2 Y],l,l); 

*/, This is ue.max LMI 

uemax_lmi = newlmi ; 
lmitermC [-uemax_lmi 1 1 Y] , 1,1); 
lmitermC [-uemax_lmi 2 1 W] , C_amp_e,M); 
lmitermC [-uemax.lmi 220], (uemax~2)); 

•/. This is ude_max LMI 

udemax_lmi = newlmi; 

lmitermC [-udemax;_lmi 1 1 Y] , eye(n) ,eye(n)) ; 
lmitermC [-udemax.lmi 2 1 W] , C_rate_e*Bbar,M) ; 
lmitermC [-udemax_lmi 2 1 Y] , C_rate_e*Abar,eye(n) ) ; 
lmitermC [-udemax_lmi 2 2 0], udemax~2) ; 

if (canardflag == 1), 

“/. This is uc_majc LMI 

ucmax_lmi = newlmi; 

lmitermC [~ucmax_lmi 1 1 Y] , 1,1); 
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lraiterm( [-ucmax_lmi 2 1 W] , C_amp_c,M); 
lmiterm( [-ucmax_lmi 220], (ucmax'2)); 



•/. This is udc_max LMI 

udcmax_lmi = newlmi ; 
lmiterm( [-udcmax_lmi 1 1 
lmiterm( [-udcmax_lmi 2 1 
lmiterm( [-udcmax_lmi 2 1 
lmiterm( [-udcmax_lmi 2 2 
end; 

lmisys = getlmis; 



Y] , eye(n) , eye(n) ) ; 

W] , C_rate_c*Bbar ,M) ; 

Y] , C_rate_c*Abar,eye(n) ) ; 
0] , udcmax~2); 



[tmin,xfeas] = f easp(lmisys , [0 200 le9 00]); 

if tmin < le-4 

'/, recover LMI variables 

y = dec2mat (lmisys, xf eas, Y) ; 
w = dec2mat (lmisys, xf eas, W) ; 

Kbar = w*M/y; 

'/, record initial condition responses 
if (canardflag == 1), 

yc = initial (Abar+Bbar*Kbar,Bbar,C_amp_c*. . . . 

Kbar , [0 0] , iO) ; 
ycmax = max(abs(yc) ) ; 

yrc = initial (Abar+Bbar*Kbar , Bbar( 1) ,... . 

C_rate_c* [Bbar*Kbar + Abar],0,i0); 

ydcmax = max(abs (yrc) ) ; 

end; 

ye = initial (Abar+Bbar*Kbar ,Bbar,C_ainp_e*Kbar, [0 0],i0); 
yemax = max(abs(ye)) ; 

yre = initial (Abar+Bbar*Kbar ,Bbar( : ,2) ,C_rate_e* ... . 
[Bbar*Kbar + Abar],0,i0); 
ydemax = max(abs(yre)) ; 
eigvalue = eig(Abar+Bbar*Kbar) ; 
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data = [data; xcg ratelimit vh split ycmax yemax .... 

ydcmax ydemax tmin utrim’ xtrim’ Kbar(l,:) Kbar(2,:) eigvalue J ] ; 



Again, for the gust penetration/state feedback problem (Version C), the code 
is very similar. Principally, the formation of the linear model (which includes gust 
dynamics) is different. After that, the code proceeds in a similar manner to the angle 
of attack/state feedback problem (Version A). That part of the code from Version C 
substantially different from Version A is shown next. 

*/, Gust penetration model formulation, 

7 . including gust dynamics and 
7 , LMI formulation. 

*/ Load appropriate data for trim 
load trimdata 

[xtrim,utrim,ytrim] = trim( ' eom_b70r ’ ,xOr ,u0 , . . . . 

[vt ; gamma] , [] , [l] , [1 ; 2] ) ; 

'/, Linear model for given cg/vh 

[af ,bf , cf ,df ] = linmod( ’ eom_b70r’ , xtrim, utrim) ; 
bf = bf (: ,2) ; df = df (: ,2) ; 



'/, Connect the actuators in series 
[A,B,C,D] = series(Aa,Ba,Ca,Da,af ,bf ,cf ,df) ; 

C= [] ; D= [] ; 

7 , Place the actuator states after the plant states 
[n,p]=size(B) ; 

T = [zeros (n-1 , 1) eye(n-l) ; eye(l) zeros (1 ,n-l)] ; 

A = T*A/T ; 

B = T*B; 

C=eye(n) ; 

D=zeros(n,p) ; 

7 , Add gust states to end of model 
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Ag=zeros(n+2 ,n+2) ; 

Ag(n+1 ,n+l)=lambdal ; 
Ag(n+2,n+2)=lambda2 ; 

Ag(l : 5,n+l)=-A(l : 5 , 2) ; 

Ag C 1 : 5 ,n+2) =-A (1 : 5 , 2) ; 

Abarl=[A zeros(n,2) ;zeros(2,n+2)] + Ag; 
Bbarl=[B;0; 0] ; 

Cbar=eye(n+2) ; 

Dbar=zeros(n+2 ,p) ; 

T = eye(n+2) ; 

T(2,n+1) = -1; 

T (2 ,n+2) = -1; 

Abar=T*Abarl/T; 

Bbar=T*Bbarl ; 



'/, Specify Static Controller Structure 
M = [eye(n) zeros(n,2); zeros (2,n+2 )] ; 

'/, Define gust i.c. 

i0=[0 000 zeros (1 ,n-4-act) zeros(l,act) -gust gust]’ 

*/, Rate/Amplitude selection matrices 

C_amp_e = 1; C_rate_e = [zeros(l,n-act) l] ; 

*/, modify for gust states 
C_rate_e = [C_rate_e 0 0] ; 

*/, A0A output matrix 

C_alpha = (57.3/vt)*[0 1 zeros(l,n-2) 0 0] ; 

V, Specify rate and amplitude constraints 

uemax = amplimit; 
udemax = ratelimit; 
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'/, form the LMIs 



setlmis ( [] ) ; 

Y = lmivar(l,[n 1;2 1]); 

W = lmivar(2,[l n+2] ) ; 

'/» This is Lyapunov lmi with real part constraint 
lyap = newlmi ; 

lmiterm( [lyap 1 1 Y] , Abar+center*eye(n+2) , 1, 's'); 
lmiterm( [lyap 1 1 W] , Bbar, M, 's'); 

'/, This is the pole placement LMI 

pplac = newlmi; 

lmiterm( [pplac 1 1 Y] , sin(theta/2)*Abar , 1, 's'); 

lmiterm( [pplac 1 1 W] , sin(theta/2)*Bbar, M, 's'); 

lmiterm( [pplac 2 1 Y], cos (theta/2) *Abar, 1); 

lmiterm( [pplac 2 1 W] , cos(theta/2)*Bbar , M) ; 

lmiterm( [pplac 2 1 -Y] , -cos(theta/2) ,Abar’) ; 
lmiterm( [pplac 2 1 -W] , -cos(theta/2)*M' ,Bbar’ ) ; 
lmiterm( [pplac 2 2 Y] , sin(theta/2)*Abar , 1, 's'); 

lmiterm( [pplac 2 2 W], sin(theta/2)*Bbar, M, 's’); 

'/, This is the Y > 0 LMI 

ypos = newlmi; 

lmiterm( [-ypos 1 1 Y],l,l); 

V. This is x_0 LMI 

xO_lmi = newlmi; 
lmiterm( [~xO_lmi 1 1 0] , 1) ; 
lmiterm( [-x0_lmi 210], iO) ; 
lmiterm( [-x0_lmi 2 2 Y],l,l); 

V, This is ue.max LMI 

uemax_lmi = newlmi ; 

lmiterm( [-uemax_lmi 1 1 Y] , 1,1); 
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lmiterm( [-uemax_lmi 2 1 W] , C_amp_e ,M) ; 
lmiterm( [-uemax_lmi 220], (uemax~2) ) ; 

'/, This is alpha_max LMI 

'/,amax_lmi = newlmi ; 

y.lmitermC [-ajnax.lmi 1 1 Y] , 1,1); 

'/,lmiterm( [-amax.lmi 2 1 Y] , C_alpha,l); 
y,lmiterm( [-amax_lmi 2 2 0], (alphalimit~2) ) ; 

'/, This is ude_max LMI 

udemax_lmi = newlmi ; 

lmiterm( [-udemax_lmi 1 1 Y] , eye(n+2) ,eye(n+2)) ; 
lmiterm( [-udemax_lmi 2 1 W] , C_rate_e*Bbar,M) ; 
lmiterm( [-udemcLX.lmi 2 1 Y] , C_rate_e*Abar ,eye(n+2) ) ; 
lmiterm( [-udemcLX.lmi 2 2 0], udemax~2) ; 



lmisys = getlmis; 

[tmin,xfeas] = feasp (lmisys , [0 200 le9 00]); 

if tmin < le-4 

y, recover LMI variables 

y = dec2mat(lmisys,xfeas,Y) ; 
w = dec2mat (lmisys ,xfeas,W) ; 

Kbar = w*M/y; 

y, record initial condition responses 

ye = initial (Abar+Bbar*Kbar , Bbar, ... . 

C_amp_e*Kbar,0,i0) ; 

yemax = max(abs(ye)) ; 

yre = initial(Abar+Bbar*Kbar,Bbar, . . . . 

C_rate_e* [Bbar*Kbar + Abar],0,i0); 

ydemax = max(abs(yre)) ; 

eigvalue = eig(Abar+Bbar*Kbar) ; 
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When the process has completed for the range of values specified for either Version 
A, B, or C, the user is returned to the MATLAB command prompt. The program 
results will have been stored in the DATA directory using the name supplied by the 
user. The user is free to view them or generate another data file. 



3. VIEWING A DATABASE 

After a database has been generated for a particular aircraft definiton, the 
program view.it. m is used to view the data. The program view.it. m presents the 
graphical user interface shown in Figure 62. 







Figure 62. Interface called by view Jt.m 

The interface is fairly self explanatory. The options are to view the three 
dimensional Tail Sizing Design Space, or to view a 2D slice of the Tail Sizing Design 
Space. Additionally, the user can specify a second database, in order to compare two 
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different aircraft definitions on a single graph. If a two dimensional view is desired, 
then the user must specify the actuator rate limit used to intersect the lower surface 
of the Tail Sizing Design Space. When ready, the user selects the button labeled 
’’Show It,” which executes the program in the GUI directory called showdata.m. 



mmmxmmmmmmxjM 



X Showdata.m 

xxxxxxxxxxxxxxx 

yi=[] ; yi_e=[] ; 



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



# / # Determine the user preferences and data file 

X names from the GUI. Go load the file from 

°/, the DATA directory. 

ipfnamel=get (loadl_ed, J string 9 ) ; 

ipf name2=get (load2_ed, J string > ) ; 

cd . . 

cd data 

eval ( [ J load * , ipf namel , 9 . mat ' ] ) ; 
cgrateA=cgrate ; 
numf iles=l ; 

if ( get (rb_f ile2, J value' ) == 1 ) , 
eval ( [ * load J , ipf name2 , 9 . mat ' ] ) ; 
cgrat eB=cgrat e ; 
numf iles=numf iles+1 ; 
end; 
cd . . 
cd model 



X This is the target actuator rate: the level plane 
# /o slicing the surface. 

udselect=str2num(get (limit_ed , * string^ ) ) ; 

X This records the user choice of either a 3D or 2D plot . 
if ( get (rb_3d, } value ') == 1) 
choice=l ; 
else , 
choice=0 ; 
end; 

X This captures the axis limits from the GUI. 
udmin=str2num(get (minrate_ed, J string ' ) ) ; 
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udmax=str2num(get(maxrate_ed, 'string')) ; 
cgmin=str2num(get (mincg_ed, ' string' ) )/100 ; 
cgmax=str2num(get(maxcg_ed, ' string' ) )/100 ; 
vhmin=str2num(get(minvh_ed, 'string')) ; 
vhmax=str2nuni(get(maxvh_ed, ’ string' ) ) ; 

*/. This captures the view point preference. 
azim=str2num(get (azim_ed, 'string' )) ; 
elev=str2num(get (elev.ed, ' string' ) ) ; 

•/. The user may wish to compare to aircraft definitions 
*/. side by side. Loop for 1 or 2 files, 
for comp=l :numf iles ; 
if (comp==l) , 
cgrate=cgrateA ; 
else, 

cgrate=cgrateB ; 
end; 

’/, Preprocess the data by parsing by tail volume. 

[m,n]=size(cgrate) ; 

num_vh=l ; 

index (num_vh)=l ; 

for i=l:m-l, 

if ( cgrate(i,l) < cgrate(i+l , 1) ), 
index (num_vh+l)=i+l ; 
num_vh=num_vh+l ; 
end; 
end; 

index (num_vh+l)=m+l ; 
f ile= ' cgrate ' ; 
for i=l:num_vh, 
p=index(i) ; 
q=index(i+l)-l; 

eval( [f ile, int2str(i) , ' = ’ ,file ,' (p:q, :);']) ; 
end; 

•/. Curve fit the rate vs. eg data for each tail 
*/. volume. The curve fit uses a simplex nonlinear 
’/, fit. The general form of the function to fit 
’/. to is two exponentials, 
global DATA 
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f ile2= ’ DATA ’ ; 
f ile3= ’ c’ ; 
lin=[] ; nlin= [] ; 
for k=l:num_vh; 

eval( [f ile2, ’ = [' ,f ile, int2str(k) , ’ ( : ,2) 

f ile , int2str(k) , ' ( : ,5)] ; ']) 

lam = [1 0] ' ; 

trace = 0; 

tol = . 1 ; 

lambda = f mins (’myf it' ,1am, [trace tol]); 

A = zeros (length(DATA( :, 1) ) , length (lambda)) ; 
for j = 1 : length (lambda) 

A(:,j) = exp(-lambda(j)*DATA( : , 1) ) ; 
end 

eval( [f ile3, ’ = A\DATA(: ,2) ; ’]) ; 
lin=[lin c] ; 
nlin=[nlin lambda]; 
end 

7, Repeat the process, this time curve fitting 
‘/, for elev amp vs. eg for each tail volume, 
global DATA_E 
file2= , DATA_E’ ; 
file3=’c’ ; 

lin_e=[]; nlin_e=[] ; 
for k=l:num_vh; 

eval( [file2, ’ = [’ ,f ile,int2str(k) , ’ ( : ,2) 

f ile, int2str(k) , ’ ( : ,6)] ; ’ ] ) 

lam = [1 0] ' ; 

trace = 0; 

tol = . 1 ; 

lambda = fmins(’myf it' ,1am, [trace tol]); 

A = zeros (length(DATA_E( :, 1) ), length(lambda)) ; 
for j = 1 : length (lambda) 

A ( : , j ) = exp (-lambda( j ) *DATA_E( : , 1) ) ; 
end 

eval( [f ile3 , ‘ = A\DATA_E( : ,2) ; »] ) ; 
lin_e= [lin_e c] ; 
nlin_e= [nlin_e lambda]; 
end 

7, Sample the functions determined above at uniform 
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7, values of eg. 
f ile4= , vh , ; 
t=cgmin: .01:cgmax; 
t=t ' ; 

for z=l:num_vh; 
for j = 1: length (t) 

r(j,z) = lin(l ,z)*exp(-nlin(l,z)*t(j)) + .... 
lin(2,z)*exp(-nlin(2,z)*t(j)) ; 
e(j,z) = lin_e(l,z)*exp(-nlin_e(l,z)*t(j)) + 
lin_e(2,z)*exp(-nlin_e(2,z)*t(j )) ; 
end ; 

eval( [f ile4, ' ( : ,z) = ’ , f ile , int2str (z) , ’ ( 1 , 1) * . . . . 
one s( length (t) , 1) ; ']) 
end; 

’/, This is the resulting mesh. It could be plotted 
'/, as is, but below we ; ll decrease the mesh size 
*/, using linear interpolation. 

ALLDATA= [] ; 

ALLDATA_E= [] ; 
for j=l:num_vh, 

ALLDATA= [ALLDATA t vh(:,j) r(:,j)]; 

ALLDATA_E= [ALLDATA.E t vh(:,j) e(:,j)]; 
end ; 

’/, Linearly interpolate the mesh between the values of 

*/, tail volume. 

ht= [] ; 

amp= [] ; 

vol=[] ; 

for k=l : ( (num_vh) -1) ; 
xi=vh(l ,k) : .01 : (vh(l ,k+l)- . 01) ; 
x=vh(l,k:k+l) ; 
for j=l : length(t) 
y=r(j,k:k+l); 
y_e=e(j ,k:k+l) ; 
yi ( : , j)=interpl(x,y,xi) ; 
yi_e( : , j)=interpl(x,y_e,xi) ; 
end 

ht= [ht yi '] ; 
amp = [amp yi_e ; ] ; 
vol= [vol ;xi '] ; 
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end 

ht=[ht r( : ,num_vh)] ; 
amp=[amp e( : ,num_vh)] ; 
vol= [vol ; vh ( 1 , num_vh) ] ; 

‘/. Plot the data, 
if (comp==l) , 
figure 
end; 

if (choice == 1) , 

'/, This is the 3D plot of the data. 
htt=ht ’ ; 

[m,n]=size(htt) ; 
for p=l:m; 
for q=l:n; 

if (htt(p,q) > udmax I htt(p,q) < 0), 
ht (q,p)=nan; 
end; 
end; 
end ; 

if (comp==l) , 
surf (t* 100 , vol ,ht ’ ) 

axis ( [cgmin* 100 cgmax*100 vhmin vhmax udmin udmax]) 

caxis([udmin udmax]) 

grid; 

view(azim , elev) 

7, This is the rate limit slice. 
ts=cgmin+ . 1 : .01 :cgmax- . 1 ; 

[mm,nn]=size(ht) ; 

pln=ones(length(ts) ,nn) *udselect ; 

C=ones(length(ts) ,nn)*2; 
surf ace (ts*100 , vol , pin’ ,C’ ) ; 
xlabel('c.g. C/»c)’) 
ylabel(’tail volume’) 
zlabel( ’ actuator rate (deg/s)’) 
else , 

surf (t* 100, vol, ht ’) 
end; 

else , 
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7 , This is the 2D intersection of the target rate limit 
7 , and surface. 
vhsc= [] ; 
cgsc=[] ; 

[m,n]=size (ALLDATA) ; 
for j=l : 3 :n 
xi = udselect; 

yi = spline (ALLDATA( : , j+2) , ALLDATA ( : , j) ,xi) ; 
vhsc=[vhsc ALLDATA(l , j+1)] ; 
cgsc=[cgsc yi] ; 
end 

7 , Fit a curve to the points representing the 2D graph. 
CGVH^cgsc’ vhsc’]; 

[p , s]=polyf it (CGVH( : , 1) ,CGVH( : ,2) ,2) ; 
lam = [1 0] ' ; 
trace = 0; 
tol = . 1 ; 

lambda = fmins( ’myf it ’ ,1am, [trace tol]); 

A = zeros (length (CGVH( 1) ), length (lambda) ) ; 
for j = 1 : length (lambda) 

A(:,j) = exp(-lambda(j)*CGVH( : , 1) ) ; 
end 

c = A\CGVH(: ,2); 

lin_s=c; 

nlin_s=lambda; 

7 , Sample the fitted function at uniform eg increments. 
vhsc= [] ; 
cgsc=[] ; 

t=cgmin : . 05 : cgmax ; 
t=t ' ; 

for j = l:length(t) 

vhsc(j,l) = lin_s (1 , l)*exp(-nlin_s (1 , 1) *t ( j) ) + 
lin_s(2, l)*exp(-nlin_s(2,l)*t(j)) ; 
end; 

cgsc=t*100 ; 

*/ Plot the 2D intersection of the level plane and surface, 
if (comp==l) , 

plot(cgsc,vhsc, ’o’ ,cgsc,vhsc) ; 
grid; 
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axis( [cgmin*100 cgmax*100 vhmin vhmax] ) ; 
xlabel('c.g. ('/.c)') 
ylabel('tail volume') 
else , 

plot (cgsc , vhsc , ,cgsc,vhsc) ; 
text(cgmax-. 15,vhmin+.025, '* = 2nd File'); 
end; 
end; 

if (numf iles==2) , 
hold ; 

clear yi yi_e r e cgrate 
end ; 

end; 

'/, Return the user to the GUI directory, 
cd . . 
cd gui 



4. THE DYNAMIC AIRCRAFT MODEL 

The model used by the Tail Sizing Design Tool simulates the longitudinal mo- 
tion of a conventional aircraft. The configuration may or may not include a canard. 
The upper shell of the dynamic model is either the Simulink file eomJ)70r.m , for 
rigid-body dynamics, or the file eomJ)70w.m , for flexible body dynamics. The trim 
and linearize routines supplied by MATLAB are use to determine the equilibrium 
point, and generate a linear model for the LMI problem formulation. The simulink 
shells above call the function eomb70r , or eomb70w , in order to compute the deriva- 
tive of the models’ states. The computation is based, in part, on the stability and 
control derivatives supplied by the user, which are retrieved via a nested function 
call explained later. The description of eornb70w follows, as the Simulink shell is self 
explanatory. 



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X function: eomb70w.m */, 

•/ 0 / •/ • / •/ •/ •/ •/ 1 / #/ •/ 0 / •/ •/ •/ •/ 1 / 1 / •/ •/ •/ •/ •/ •/ •/ • / 0 / •/ 0 / 0 / 



165 



% The vector passed to the function is the 
% vehicles' states, 
function accel = eomb70w(x) 

% Retrieve the center of gravity , tail volume 
7. and canard volume parameters 

load cgvh; 

% Retrieve the normalized aircraft derivatives 
•/. at the center of gravity and control volumes 
i, specified. Retrieve the flexible mode frequency. 

[rho , a , CDO , CL , CM , QF , omegai ,mi,Ti,s,b,c,m,To, Iyy , lh , 1 c] 
b70dataw(xcg,vh, split) ; 
load flexdata 

it Parse the vector passed to the function into 
it its components and use them to compute terms 
i, necessary to "denormalize" the force and moment 
it coefficients. 

dc = x(l); de = x(2); dtrt = x(3); 

u = x(4) ; w = x(5) ; q = x(6); 

theta = x(7); eta = x(8); etad = x(9); 

it Calculate the total velocity, vt , and dynamic 
it pressure gravity (assumed constant) and 
it angle of attack. 



vt 


= 


(u~2 + w~2) ~ . 5 


qbar 


= 


. 5*rho* (vt~2) ; 


g 


= 


32.2; 


alpha 


= 


asin(w/vt) ; 


norml 


= 


rho*vt~2*s/2 ; 


norm2 


= 


rho*vt*s*c/4; 


norm3 


= 


rho*vt~2*s*c/2 


norm4 


= 


rho*vt*s*c‘2/4 


norm5 


= 


rho*vt“2*s*c/2 


norm6 


= 


rho*vt*s*c~2/4 



it Based on the stability and control derivatives 
it retrieved by the function {\em b70dataw.m>, the 
it coefficients above and the vehicles' states, 
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•/. the total lift, drag, pitching moment and elastic 
% generalized force on the vehicle is computed. 

L = norml* (CL ( 1 ) + CL(2)*alpha + CL(7)*de + CL(8)*dc + 

CL(5) *eta) + norm2* (CL(3) *0 + CL(4)*q + CL(6)*etad); 

Ln = L/ (qbar*s) ; 

Dn = (CDO + 1.01*Ln~2/pi/a); 

D = Dn*qbar*s; 

M = norm3* (CM( 1) + CM(2)*alpha + CM(7)*de + CM(8)*dc + 

CM(5)*eta) + norm4* (CM(3)*0 + CM(4)*q + CM(6)*etad); 

Qeta = norm5*(QF(l) + QF(2)*alpha + QF(6)*de + QF(7)*dc .... 

+ QF(4)*eta) + norm6* (QF(3)*q + QF(5)*etad); 

•/o The lift, drag and gravity force are resolved in 
•/. the body reference frame. Then, the linear, angular 
•/. amd elastc accelerations are computed in the body 
°/ reference frame. 

udot = -q*w - g*sin(theta) + (-cos (alpha) *D 
sin(alpha)*L + To*dtrt)/m; 

wdot = q*u + g*cos(theta) + (sin(alpha)*D - .... 
cos (alpha) *L)/m; 
qdot = M/Iyy; 

etadd = -omegai(l) ~2*eta - 2*omegai*Ti*etad + Qeta/mi; 

Notice that the prior program utilized a nested function call to retrieve the 
total vehicle stability and control derivatives. These derivates are supplied by the 
function blOdataw , which provide user supplied values of the wing-body, tail and 
canard contributions to lift and pitch. If a flexible mode is retained, the user must 
also specify the mode shape, in vacuo frequency, and generalized mass of the first 
bending mode. Presented below is data for the generic HSCT model termed Ref A. 
When the flexible mode is included, the model is termed Ref B. 



«/ o/ •/ •/ •/ •/ •/ •/ •/ •/•/ •/ •/•/ •/•/ •/ •/ «/ •/ •/•/ •/ •/ •/ v •/ •/ 



• f% u ft tt ft ft ft ft ft ft ft ft it it it u u u ft ft ft ft ft ft tt ft ft 



'/, function: b70dataw.m '/, 
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•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/ •/ '/ V •/ V •/ •/ V •/ V •/ •/ •/ e 



'«/•/•/«/•/•/•/• 



‘/, To begin, the user roust specify the physical 
*/, description of the aircraft including the 
'/, frequency and generalized mass associated 
'/, with the flexible mode. 



'/, parameter/value/units/description 
G = 32.174; '/, gravitational constant 



g = 


G; ’/. fps~2 


w = 


480425; */. lbs 


m = 


W/G; */. slugs 


iyy = 


= 2 . 181e7 ; */. slug-ft~2 


s = 


6298; •/. ft '2 


c = 


78.5; •/. ft 


cbar 


= c; '/. c 


b = 


105; */. wing span 


U = 


422; */. fps 


rho = 


= .0023769; '/, density 


Q = 


.5*rho*U*U; '/, lbf/ft*ft 


a = 


1.751; % aspect ratio 


lh = 


85; '/, tail moment arm ■ 



lc = 85; ’/, canard moment arm to ac 

ih = 0; ‘/, tail incidence angle 

To = 1000; 7 , nominal thrust lbs 

cloc = .05; */, canard location 

tloc = .95; '/, aero center of tail location 

vc = split; */, canard volume 

vh = vh; '/, tail volume 

omegai = 6.29; '/. flexible frequency 

mi = [500] ; '/, generalized mass 

Ti = .02; */, flexible damping 

% The stability and control derivatives of the 
7. canard, wing-body and tail are entered. Those listed 
•/. below are for a generic HSCT "like" aircraft termed 
*/, ref A. Following the development of Wykes 'method, 
the user can enter terms that residualize the effects 
% of the elastic motion. In this case, only the rigid 
’/, body derivatives have been entered. The nomenclature 
% uses either no extension or the extension (_r) to 
'/, indicate that the term is a rigid body term. The 
# /, extension (_e) is used for flexible body terms. 
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*/, Finally, the extension (_eor) is used for terms that 
*/, scale the rigid body terms due to flexible effects. 



'/, Drag data 
CDO = 0.002; 

CDU = 0.0; 

'/, Lift data 
clift_total = 0.161; 
clift_vh_term = 0.801642; 
clift_wb_af 0_r = 0.04; 
clift_dwb_af 0_er = 0.0; 
clift_wb_af _r = 0.0422; 
clift_wb_af _eor = 1.0; 
clift_wb_q_r = 0.0558; 
clift _wb_q_eor = 1.0; 
clif t_wb_adf _r = 0.0; 
clift _wb_adf _eor = 0.0; 
clift_wb_n_e = 0.0; 
clift_wb_dq_e = 0.0; 
clift_h_q_r = 0.123; 
clift_h_q_eor = 1.0; 
clift_h_adf_e = 0.0; 
clift_h_n_e = 0.0; 
clift _h_dq_e = 0.0; 
clif t_h_ah0_e = 0.0; 
clift_h_ah_r = 0.0455; 
clift_h_ah_eor = 1.0; 
clift _e_af0_r = 1.8; 
clift_e_af 0_er = 0.0; 
clift_dah_dloadh = 0.0; 
clift_e_af_r = 0.71; 
clift_e_af _eor = 1.0; 

*/, Lift Data - Canard 
clift_c_q_r = 0.123; 
clift_c_q_eor = 1.0; 
clift_c_adf _e = 0.0; 
clift_c_n_e = 0.0; 
clift_c_dq_e = 0.0; 
clif t_c_ac0_e = 0.0; 
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clift_c_ac_r = 0.0455; 
clift_c_ac_eor = 1.0; 
clift_dac_dloadc = 0.0; 

CLU = 0; 

'/, Pitching moment data 
cm_total = -0.063618; 
cmO_wb_r = 0.0; 
cmO_wb_er - 0.0; 
dcm_dcl_wb_r = 0.02; 
dcm_dcl_wb_er = 0.0; 
cm_wb_q_r = -0.0298; 
cm_wb_q_eor = 1.0; 
cm_wb_adf_r = -0.0067; 
cm_wb_adf _eor = 1.0; 
cm_wb_n_e = -0.01; 

cm_wb_dq_e = 0.0; 

CMU =0.0; 

'/, Based on theses derivatives and the physical 
'/, characteristics of the aircraft, some useful 
'/, coefficients are calculated. 

kl = cbar/(2*U) ; 
k2 = cbar*vh/lh; 
k3 = S*cbar/lh; 
k2_c = cbar*vc/lc; 
k3_c = S*cbar/lc; 

'/o Tail contribution coefficients 

Alpha. 1 = clift_h_ah0_e. . . 

/ (l-(clift_h_ah_r*clift_h_ah_eor*clift_dah_dloadh*Q*vh*k3) ) ; 
Alpha_2 = -(clift_h_ah_r*clift_h_ah_eor*. . . . 

(clift_e_af 0_r+clift_e_af 0_er) ) . . . 

/(l- (clift_h_ah_r*clift_h_ah_eor*clift_dah_dloadh*Q*vh*k3) ) ; 
Alpha.O = Alpha. 1 + Alpha_2; 

Alpha_ih= clift _h_ah_r*clift_h_ah_eor. . . 

/ (l-(clift_h_ah_r*clift_h_ah_eor*clift_dah_dloadh*Q*vh*k3)) ; 
Beta = (clift_h_ah_r*clift_h_ah_eor* . . . . 

(l-clift_e_af _r*clift_e_af_eor) ) . . . 

/ (1- (cl if t_h_ah_r*clift_h_ah_eor*clift_dah_dloadh*Q*vh*k3) ) ; 
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°/ 0 Canard contribution coefficients. Note, the 
*/ 9 canard is treated simply as a forward tail. 

Alpha_0c= clif t_c_acO_e . . . 

/(l-(clift_c_ac_r*clift_c_ac_eor*clift_dac_dloadc . .. . 
*Q*vc*k3_c) ) ; 

Alpha_ic= clif t_c_ac_r*clif t_c_ac_eor . . . 
/(l-(clift_c_ac_r*clift_c_ac_eor*clift_dac_dloadc. . . . 
*Q*vc*k3_c) ) ; 

Beta_c = (clift_c_ac_r*clift_c_ac_eor*(l-clift_e_af _r . . . . 
*clift_e_af _eor) ) . . . 

/ (1- (clif t_c_ac_r*clift_c_ac_eor*clift_dac_dloadc*Q*vc*k3_c) ) ; 

•/. Finally, the total vehicle stability and control 
*/, derivatives are computed based on the canard, 
y. wing-body and tail contributions. The build up 
# /. is based on a standard development presented in 
7 . many aircraft dynamics texts; here I used Etkin. 

°/, The computation of the elastic stability and control 
# /, derivatives follow the development earlier in this text. 

# /o Form static stability derivaties needed for trim. 

CLO = clif t_wb_af 0_r + clif t_dwb_af 0_er + .... 
k2*Alpha_0 + k2_c*Alpha_0c; 

CLA = 57.3*(clift_wb_af_r*clift_wb_af_eor + .... 
k2*Beta + k2_c*Beta_c) ; 

CLIS = k2*Alpha_ih; 

CLIC = k2_c*Alpha_ic ; 

CMO = cmO_wb_r + cmO_wb_er + (dcm_dcl_wb_r + .... 
dcm_dcl_wb_er) *(clif t_wb_af 0_r + clif t_dwb_af 0_er) .... 

- Alpha_0*vh +Alpha_Oc*vc + CL0*(xcg- . 25) ; 

CMA = 57.3*((dcm_dcl_wb_r + 

dcm_dcl_wb_er)*(clift_wb_af _r*clif t_wb_af .eor) . . . 

- Beta*vh + Beta_c*vc) + CLA*(xcg- . 25) ; 

CMIS = -Alpha_ih*vh + CLIS* (xcg- . 25) ; 

CMIC = Alpha_ic*vc + CLIC* (xcg- . 25) ; 

# / 8 Compute dynamic stability derivatives needed 
°/, for linearization. 
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CLQ = 57 . 3* (clif t_wb_q_r*clift_wb_q_eor + ... 
clif t_h_q_r*clif t_h_q_eor*k2 - clif t_c_q_r* . . . . 
clift_c_q_eor*k2_c) ; 

CLAD = 57 .3*kl*(clif t_wb_adf _r*clift_wb_adf_eor + ... 
clift_h_adf _e*k2 + clif t_c_adf _e*k2_c) ; 

CMQ = 57.3*((dcm_dcl_wb_r + dcm_dcl_wb_er)* . . . . 
(clift_wb_q_r*clift_wb_q_eor) + cm_wb_q_r*cm_wb_q_eor -.... 
clif t_h_q_r*clif t_h_q_eor*vh 

clif t_c_q_r*clif t_c_q_eor*vc) +CLQ* (xcg- . 25 ) ; 

CMAD = 57 . 3* ( (dcm_dcl_wb_r + dcm_dcl_wb_er)* . . . . 
(cm_wb_adf_r*cm_wb_adf _eor) . . . 

- clift_h_adf _e*vh - clift_c_adf _e*vc) + CLAD* (xcg- . 25) ; 

"/, Compute stability derivatives for the 
’/. generalized elastic coordinates. When these are 
Y, non-zero, the model is aeroelastic. The resulting 
’I, model is termed Ref B is Chapter 2. 

[phi , phidot] =modeshpe ( . 25*250) ; 

CLeta = CLA*phidot; 

CLetad = CLA*phi/U; 

CMeta = CLA* (xcg- .25)*phidot ; 

CMetad = CLAD*phi*(xcg- .25)/U; 

Q0 = -CL0*phi; 

QA = -CLA*phi; 

QQ - -CLQ*phi; 

QETA = -CLA*phidot*phi ; 

QETAD = -CLAD*phi*phi/U; 

[phi , phidot] =modeshpe (tloc*250 ) ; 

QIS = -CLIS*phi; 

[phi .phidot] =modeshpe(cloc*250) ; 

QIC = -CLIC*phi ; 

*/, Combine the individual terms into a vector 
y. to be returned. 

CL = [CLO CLA CLAD CLQ CLeta CLetad CLIS CLIC] ; 

CM = [CMO CMA CMAD CMQ CMeta CMetad CMIS CMIC] ; 

QF = [Q0 QA QQ QETA QETAD QIS QIC] ; 
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APPENDIX B. RFTPS SUPPORT 



1. INTRODUCTION 

This appendix documents the C-code unique to the RFTPS design that is not 
automatically generated via the autocode tools. The focus is on the serial module, 
which processes both the IMU and GPS data. More complete hardware descriptions 
of the I/O modules supported by the AC100 architecture can be found in the AC100 
Model C30 Supplemental Reference Manual. Wiring diagrams and pin assignments 
are on file in the Avionics Lab. 



2. SERIAL MODULE 

The serial module contains two multi-mode, bi-directional, channels. The user 
interface to the device is through the file userser.c. Within user.ser.c , three functions 
must be defined. The first, get_SERIAL_parameters, is used to initialize the channels. 
The second, user_SERIAL_out, defines the output from each channel and the third, 
user_SERIALJn, processes received data in the input buffer. Currently, the serial 
module is not used to send data. Some of the shell in userser.c is not shown. The 
documentation focuses on RFTPS specific code associated with processing the GPS 
and IMU data. 

/ *************************** ********************************* 

** File : user.ser.c 

** Project : RFTPS 
** Edit level : 4 
** 

** Abstract: : File contains functions which the user 

** must define to interface with IP-SERIAL device 
** driver. 

** The functions must be called 

** get_SERIAL_parameters 
** user_sample_SERIAL_in 
** user_SERIAL_out 
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** 

** Templates for the functions are provided. 

** 

** get_SERIAL_parameters 

** Function sets the asynchronous communication 
** parameters for the IP-SERIAL module. Ring buffer 
** sizes used to store received data must also be 
** specified. 

** 

** user_SERIAL_out 

** Function is called every scheduler interval. The 
** user is responsible for creating a byte stream from the 
** models floating point outputs. The user must ensure 
** that the when writing these bytes to the output buffers 
** that the buffers are not overflown. 

** 

** user_sample_SERIAL_in 

** Function is called every sampling interval. The 
** user is responsible for filling the floating-point 
** vector which is used as input to the model for 
** the current sampling interval. 

** 

** Modifications: 

** 



** 


Creation 


7-1—93 


Henry Tominaga 






** 


Revised 


8-23-93 


Brent Roman 






** 


Revised 


11-18-93 


Steve Lynch 






** 


Revised 


9-1-94 


Steve Lynch 






** 


Revised 


1-17-96 


Eric Hallberg 


[New 


IMU] 


** 


Revised 


3-15-96 


Eric Hallberg 


[IMU 


Serial A 


** 






/ 


GPS 


Serial B] 


** 


Revised 


: 5-01-96 


Eric Hallberg 


[IMU 


binary] 



*/ 



#include 

#include 

Jfinclude 

#include 

#include 

#include 

#include 



<stddef .h> 
<stdlib.h> 
"sa_types .h" 
"stdtypes .h" 
"actypes .h" 
"ioerrs .h" 
"errcodes .h" 
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#include "iodefs.h" 

#include "iodriver.h" 

#include "ipserial.h" 

#include "printx.h" 

#def ine NULL 0 

/* semaphores and serial parameters for each 
physical channel */ 
private struct 

{ ISI.BOOLEAN allSent ; ISI.BOOLEAN broken; 
unsigned baud;} line_state [2] ; 

struct user_type 

int update_interval ; 
int update_count ; 



>; 



typedef struct _bytes 

unsigned bytel :8; 
unsigned byte2 :8; 
unsigned byte3 :8; 
unsigned byte4 :8; 

}_bytes ; 

typedef union float_char 

float fl; 

_bytes ch; 

} float_char; 

/* global variables used in user_ser_in ch A 
IMU data processing */ 

u_int buff er_data[40] ; 

int index ; 

int tick = 0; 

int f irst_frame_a = 0; 
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float last_f loat_a[l4] ; 
float last_f loat_a_l [14] ; 
float last_f loat_a_2 [14] ; 
float last_f loat_a_3 [14] ; 
int missed_cr = 0; 



/* global variables used in user_ser_in ch B 
GPS processing */ 

int f irst_f rame_b = 0; 

float last_f loat_b[50] , sol=299792458.0, Ilf 0=1575420000 . 0 , k2; 
unsigned long icpoldl, icpold2, icpold3, icpold4, 
icpold5, icpold6; 
int num_bytes_old=0 ; 

/* Function: get_SERIAL_parameters +++++++++++++++++++++ 

** Returns : 

** NONE 
*/ 

public void get_SERIAL_parameters 

( 

unsigned int hardware_channel , 

volatile struct user_param *device_param, 

volatile struct ring_buf f er_param *rec_buff er , 

IOdevice ^device) 

{ 

struct user_type *user_ptr; 
serial_param_type *serptr; 

serptr = device->parameters; 

if (SERIAL_USER_PTR == NULL) 

{ 

SERIAL_USER_PTR = (void *)malloc (sizeof (struct user_type)); 
user.ptr = (struct user_type *)SERIAL_USER_PTR; 
user_ptr->update_interval = 1000; 
user_ptr->update_count = 0; 

} 



/* Both IMU and GPS sensors send their data at 



176 



9600 baud/RS-232 standard */ 

if (hardware_channel == chanA | | hardware_channel == chanB) 

{ 

/* set parameters for channel A (IP-SERIAL channel 1)*/ 
device_param->parity = NONE; 

device_param->baud_rate = 9600; 

device_param->stop_bits = ONE; 

device_param->transmit_data_size = 8; 
device_param->receive_data_size = 8 ; 
device_param->clock_multiplier = 16 ; 

/* Set size for receive ring buffer. 

Here it is set to lOx the message size */ 
rec_buf f er->buff er_size = 600; 

> 

else 

{ 

printxC" INVALID CHANNEL\n") ; 

} 

} /* get_SERIAL_parameters */ 



/* Function: user_SERIAL_out 
*/ 



public RetCode user_SERIAL_out (IOdevice ^device, 
float model_f loat [] , 
u_int ser_channel) 



Byte 

u_int 

f loat_char 
RetCode 



cbuff er [20] ; 

i; 

data; 

return_val; 



return_val = OK; 

/*************************** it:***********!*:**)*:**)*: 

* Given floating point model output, 

* create buffer which contains bytes 

* to be transmitted across 

* serial channel. 

if (numbytes_in_buff er (device->parameters) == 0) 

{ 

for (i = 0; i < 5; i++) 
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{ 

data.fl = model_f loat [i] ; 
cbuffer[0] =(Byte) data . ch . bytel ; 
cbuffer[l] =(Byte) data. ch . byte2 ; 
cbuffer[2] =(Byte) data. ch.byte3; 
cbuffer[3] =(Byte) data.ch.byte4; 

* Fill the output buffer with data to be transmitted 

* by the background portion of the serial driver 

return_val = write_serial(device->parameters ,cbuff er ,4) ; 
if (return_val == -1) 

{ 

printx("Error : Buffer overflow in User_ser on 
output to serial\n") ; 

> 

> 

> 

return OK; 

} /* user_SERIAL_out */ 

/* Function: user_sample_SERIAL_in +++++++++++++++++ 

*/ 

public RetCode user_sample_SERIAL_in(IOdevice ^device, 
float model_f loat [] , 

u_int ser_channel) 

{ 

/*********!tcREADING CHANNEL A ****** IMU ********/ 
if (ser_channel==chanA) 

{ 

/* Create a software buffer to hold receive 
buffer data */ 
u_int soft_buff er [600] ; 

/* Declare varibiables needed */ 
unsigned char item[l] ; 
float k; 
int y ,z, j ; 
int s = 0; 
int p = 0; 
u_int w_high; 
u_int w_low; 
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struct user.type *user_ptr; 

/* The IMU data is sent as twos complement 
and must be scaled*/ 

float scale [] = {180/8191.0, -2.0/8191.0, 

-2.0/8191.0, -2.0/8191.0, 
200.0/8191.0, 200.0/8191.0, 

200.0/8191.0, 180.0/8191, 
180.0/8191.0, 10.0/8191.0, 

10.0/8191, 10.0/8191.0, 

10.0/8191.0, 10. 0/8191. 0>, 

def ault_data[] = { 10.0, 10.0, -10.0, 10.0, 10.0, 10.0, 
10.0, 10.0, 10.0, 1.0, 1.0, 1.0, 1.0, 1.0}; 

/ ********** ********* * ************************** 

* set user pointer to buffer allocated 

* by get parameters 

* this buffer is passed around with the 

* structure device 

* and should only be accessed via the SERIAL_USER_PTR 

* 

*********************************************/ 
user_ptr = SERIAL_USER_PTR; 

/* The first time the serial port is read; 

load default data */ 
if (f irst_f rame_a==0){ 
for (j=0; j<14; j++){ 
last_f loat_a[j] = def ault_data[j] ; 
last_f loat_a_l [j] = def ault_data[j] ; 
last_f loat_a_2 [j] - def ault_data[j] ; 
last_f loat_a_3 [j] = def ault_data[j] ; 

}/* end if*/ 

model_float = last_f loat_a; 
index = 0; 
f irst_f rame_a = 1; 
return OK; 

}/* end if*/ 

/* Process hardware buffer by moving data into 



179 



the software buffer */ 

while ( (numbytes_in_buff er(device->parameters) ) >1 ){ 

read_serial (device->parameters , 1 , item) ; 
sof t_buf f er [s] =(u_int) item[0]; 
s=s+l ; 

} /* end while */ 

/* The IMU data string ends with a carriage return. 

Always sync on the carriage return*/ 
while ( p < (s)){ 
if (soft_buff er [p] == '\r'){ 
missed_cr = 1 ; 
index = 0; 

}/* end if*/ 
else-C 

if (missed_cr == l){ 

buffer_data[ index ] = soft_buffer [p] ; 

}/* end if */ 

index=index + 1; 
if (index == 30) { 
index = 0; 
missed_cr = 0; 

}/*end if*/ 

}/* end else*/ 
p=p+l; 

} /* end while */ 

/* Implement some rudimentary error detection to remove 

* blatantly bad data from the IMU. Since the IMU does not 

* send a check sum, use multiple buffers to compare present 

* with past values . Remove data readings that are not 

* physically possible. Send last good data. */ 

tick = tick + 1; 
if (tick == 4) 
tick = 1; 



if (tick == 1) 

{ 

for ( y=0 ; y<14; y++){ 
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/* The number is in binary, 2's complement. 

Convert to decimal and scale */ 
w_high = buff er_data[y*2] ; 

w_low = buf f er_data[y*2+l] ; 
k = 128*(w_high & 0x7F) + (w_low & 0x7F) ; 

if (k > 8191) 
k = k - 16384; 

last_f loat_a_l [y] = k*scale[y]; 

if (abs (last_f loat_a_l [y] - last_f loat_a_3 [y] ) < .1) 
last_f loat_a[y] = last_f loat_a_l [y] ; 

} /* end for loop */ 

}/*end if*/ 

if (tick == 2) 

{ 

for ( y=0 ; y<14; y++){ 

w_high = buff er_data[y*2] ; 

w_low = buff er_data[y*2+l] ; 
k = 128* (w_high & 0x7F) + (w_low & 0x7F) ; 

if (k > 8191) 
k = k - 16384; 

last_f loat_a_2 [y] = k*scale[y]; 

if (abs (last _float_a_2 [y] - last_f loat_a_l [y] ) < . 1) 
last_float_a[y] = last_f loat_a_2 [y] ; 

} /* end for loop */ 

}/*end if*/ 

if (tick == 3) 

{ 

for ( y=0; y<14; y++){ 

w_high = buff er_data[y*2] ; 
w_low = buff er_data[y*2+l] ; 
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k = 128*(w_high & 0x7F) + (w_low & 0x7F) ; 

if (k > 8191) 
k = k - 16384; 

last_f loat_a_3 [y] = k*scale[y] ; 

if (abs(last_float_a_3[y] - last_f loat_a_2 [y] ) < .1) 
last_f loat_a[y] = last_float_a_3[y] ; 

} /* end for loop */ 

}/*end if*/ 

/* Pass the IMU data back up. 

for (j=0; j<14; j++){ 

model_f loat [j] = last_float_a[j] ; 

}/* end for*/ 

}/* end if ch a */ 

/****** READING CHANNEL B**** GPS ********/ 
if (ser_channel==chanB){ 

unsigned char itemb[l], msgl [600] , checksuml [l] , 
checksum2 [1] , my checksuml, my checksum2, status, 
model, mode2, mode3, mode4, mode5, mode6, 
msg3[45], checksum3 [1] , mychecksum3, msg4[62] , 
checksum4[l] , mychecksum4; 

struct user_type *user_ptr; 

unsigned long i, j, icpl, icp2, icp3, 
icp4, icp5, icp6, 
ephl, eph2, eph3, eph4, 
eph5 , eph6 ; 

long lat , Ion, vel, heigps, heimsl; 

int numl, num2, num3, satl, sat2, 
sat3, sat4, sat5, sat6, 
sattl, satt2, satt3, satt4, 
satt5, satt6, num4, 
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velnor, velup, veleast; 



float psrngl, psrng2, psrng3, psrng4, 
psrng5, psrng6, psrngratel, psrngrate2, 
psrngrate3, psrngrate4, psrngrate5, 
psrngrate6 , loctime , sattimel , satt ime2 , 
sattime3, sattime4, sattime5, sattime6, 
ecefx, ecefy, ecefz, velnorl, veleastl, 
velupl, hd, timeref, pscorl, pscor2, 
pscor3, pscor4, pscor5, pscor6, 
psrcorl, psrcor2, psrcor3, psrcor4, 
psrcorS, psrcor6, locsec, locsecf, 
satsecl, satsec2, satsec3, satsec4, 
satsec5 , satsec6 , sat secf 1 , satsecf 2 , 
satsecf3, satsecf 4, satsecf 5, 
satsecf6, heigpsl, heimsll, 
veil, latl, lonl, ecefxl, ecefyl, ecefzl; 

/********************************************************** 

* set user pointer to buffer allocated by get parameters * 

* this buffer is passed around with the structure device * 

* and should only be accessed via the SERIAL_USER_PTR * 

* define * 

user.ptr = SERIAL.USER.PTR; 

* set user pointer to buffer allocated by get parameters * 

* this buffer is passed around with the structure device * 

* and should only be accessed via the SERIAL. USER. PTR * 

* define * 

**********************************************************/ 

/* If first time read, load default data */ 
if (f irst.f rame_b==0) { 
f or (i=0 ; i<50 ; i++) { 
last.f loat.b [i] =7.0; 
f irst .frame. b=l ; 

icpoldl = 0 . 0 ; icpold2=0 .0; icpold3=0 . 0 ; icpold4=0 . 0 ; 
icpold5=0 . 0 ; icpold6=0 . 0; 

}/*end for*/ 
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}/*end if*/ 



/* The message is 61 bytes long and comes in once a second. 
* Wait for the whole message to come in. */ 
numl = (numbytes_in_buffer(device->parameters) ) ; 
if ((numl > 60) k (numl == num_bytes_old) ){ 

for (i=0; i<numl ; i++){ 

read_serial (device->parameters , 1 , itemb) ; 
msgl [i]=itemb[0] ; 

}/* end for */ 

}/*end if message in*/ 

/* The message is in ASCII Hex and starts with ©@BA. 

Search for start of string. */ 
for (i=4; i<numl ; i++){ 

if (msgl [i-4]= =, @' k msglCi-S] 3 ^®' & 
msgl [i-2]=='B' & msgl[i-l]== , a' ){ 

checksuml[0] = msgl [i+61] ; 

/♦Convert message bytes to decimal numbers*/ 

lat = msgl [i+11] *0x1000000 + msgl [i+12] *0x10000 + 
msgl [i+13] *0x100 + msgl[i+14]; 

Ion = msgl [i+15] *0x1000000 + msgl [i+16] *0x10000 + 
msgl [i+17] *0x100 + msgl[i+18]; 

heigps =msgl [i+19] *0x1000000 +msgl [i+20]*0xl0000 + 
msgl [i+21] *0x100 + msgl [i+22] ; 

heimsl = msgl [i+23] *0x1000000 +msgl [i+24] *0x10000 + 
msgl [i+25] *0x100 +msgl[i+26]; 

vel = msgl [i+27] *0x100 + msgl[i+28]; 

hd = msgl [i+29] *0x100 + msgl [30]; 

hd = hd/10.0; 
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status = msgl[i+60]; 

vell=vel/lOO . 0 ; heigpsl=heigps/lOO . 0 ; heimsll=heimsl/100 . 0 ; 
latl=lat/100 . 0 ; lonl=lon/lOO . 0 ; 

/*GPS uses a check sum. Check if message is correct*/ 
mychecksuml= ’ B ’ " ’ a 1 ; 
for(j=0; j<61;j++){ 



mychecksural~=msgl [i+j] ; 

}/*end for*/ 

/* If its correct, update GPS data, 
if (mychecksuml == checksum!. [0] ){ 

last_f loat_b [0] =latl ; 
last_f loat_b [l] =lonl ; 
last_f loat_b [2] =heigpsl ; 
last_f loat_b [3] =heimsll ; 
last_f loat_b [4] =vell ; 
last_f loat_b [5] =hd; 
last_f loat_b [6] =status ; 
last_f loat_b [7] =checksuml [0] ; 
last_f loat_b [8] =mychecksuml ; 

}/*end if checksum ok*/ 

i=i+60 ; 

}/* end if @@Ba */ 

/* Pass the GPS data back up */ 
for (j=0; j<8; j++){ 
model_f loat [j] = last_float_a[j] ; 
>/* end for*/ 

}/* end if ch b */ 
return OK; 

} /* user_sample_SERIAL_in */ 
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